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EXECUTIVE SUMMARY 


Introduction 


The bulk-flow analysis results for this contract are incorporated in the following publications: 

.(a) Childs, D. (1991a), "Fluid-Structure Interaction Forces at Pump-Impeller- 
^ Shroud Surfaces for Axial Vibration Analysis," ASME Trans., Journal of 

Vibration and Acoustics, Vol. 113, pp. 113-115, January 1991 

(bf Childs, D. (1991b), "Centrifugal Acceleration Modes for Incompressible Fluid 
^ in the Leakage Annulus Between a Shrouded Pump Impeller and Its Housing," 

ASME Trans., Journal of Vibration and Acoustics, Vol. 1 13, pp. 209-218, 

April 1991 

y^/3 (c) JVilliams, J. and Childs, D. (1992aX "Influence of Impeller Shroud Forces on 
Pump Rotordynamics," ASME Trans., Journal of Vibration and Acoustics, Vol. 
113, pp. 509-515, October 1991 

(df Childs, D. (1992b), "Pressure Oscillation in the Leakage Annulus Between a 
U Shrouded Impeller and Its Housing Due to Impeller-Discharge-Pressure 

Disturbances, ASME Trans., Journal of Fluids Engineering, V ol. 114, pp. 61- 
67, March 1992 

( e ) Cao, N. (1993), "Compressibility Effects on Rotor Forces in the Leakage Path 
Between a Shrouded Pump Impeller and Its Housing," M.S. M.E. Thesis, Texas 
A&M University, August 1993 

Computational Fluid Mechanics (CFD) results developed by Dr. Erian Baskharone are 
reported separately. Two copies of these publications, were submitted as preliminary final 
reports under the terms of this contract and are incorporated in this final report. 

The results of this study and the publications above can be summarized as follows: 


Impeller Forces for Axial Vibration Analysis 


Initial bulk-flow analyses for impellers considered radial reaction forces developed by 
impellers due to lateral rotor motion Reference (a) above examined the axial reaction forces, 
concluding that no resonance exists for axial motion. Methods are presented for calculating 
stiffness, damping, and added mass coefficients for axial vibrations of turbopumps. 


Centrifugal Acceleration Modes 

Prior to this study, Childs (1989) calculated reaction force components for an impeller due to 
precession of a pump rotor at nonsynchronous frequencies. His results showed unexpected 
peaks in the force components which he ascribed to fluid "resonance", arising from the 
centrifugal-acceleration term in the momentum equation. Dr. Brennen at Cal Tech questioned 
this terminology, suggesting that the term resonance could only be supported if an analysis 
confirmed that the governing system actually had complex eigenvalues at or near the locations 
of peak amplitudes. Childs (1991b) cited above, yielded roots and eigenvectors at the 
predicted locations, confirming the prior predictions of fluid resonances. 


Influence of Impeller Shroud Forces on Pump Rotordynamics 

A question presented by the initial predictions that impeller reaction forces could contain 
"peaks" was: How should radial and circumferential reaction forces be modeled if they can 
not be modeled with stiffness, damping, and added-mass coefficients? Furthermore, what 
influence do the predicted peaks have on rotordynamics? Williams and Childs (1992a) 
developed linear and nonlinear analysis procedures for incorporating the frequency-dependent 
radial and circumferential force coefficients into a rotordynamic analysis. Transient nonlinear 
analysis used the predicted reaction forces directly as a function of the instantaneous 
normalized precession frequency. Synchronous response due to imbalance proceeds directly, 
replacing the radial and circumferential force coefficients with direct and cross-coupled 
stiffness coefficients which are a function of running speed. Eigenanalysis is iterative, since 
the stiffness coefficients depend on the precession frequency which is in turn defined by the 
eigenvalue. 

The predicted impeller-force peaks were shown to have a major influence on rotordynamics 
for the model considered. However, the results are difficult to generalize to turbomachinery 
rotordynamics. 


Pressure Oscillation Excitation of the Flow in an Impeller Shrouded Annulus 

The analyses cited above all considered reaction forces due to impeller shroud motion. 

Childs (1992b) looked at the separate problem of pressure and velocity predictions in the 
annulus due to precessing pressure oscillations at the discharge of the impeller. This analysis 
was not aimed at rotordynamics . The impeller is assumed to be spinning but not precessing. 


Pressure oscillation from the impeller discharge propagate down through the annulus and out 
the exit wearing-ring seal. The analysis showed that multiple sharp resonances could be 
excited in the annulus and that they typically resulted in peak pressure oscillations at or near 
the exit wearing-ring seal. These results suggested a possible explanation for the internal 
cooling problem for the KEL-F exit seals of the SSME-HPFTP. However, the applicability of 
the analysis was limited because an incompressible-flow model was used. 


Compressibility Effects on Rotordynamics and Leakage and Pressure in an Impeller 
Annulus 


Nhai’s (1993) thesis extended previous models by incorporating fluid compressibi lty. ai 
uses a barotropic model for which the viscosity and density are functions of the pressure 
(only) Adding compressibility to the model means that acoustic modes can be generated in 
addition to the "centrifugal-accelerations" modes which were present 

Acoustic analysis normally discards perturbation terms which are included in Nhai s general 
perturbation analysis. 

Nhai used the HPFTP first stage impeller for his analysis. He analyzed the exit weanng-nng 
seal leakage-AP relationship using a code developed by Momson et al. (1983). Fluid 
properties were modeled via an NBS code, McCarthy et al. (1986). 

Nhai examined the influence of compressibility on both rotordynamic characteristics of pumps 
and pressure and flow oscillations within the annulus due to pressure perturbations of the 
impeUer exit. Concerning rotordynamic-response characteristics, Nhai s analysis sought to 
answer the following basic questions! 


(a) What influence does compressibility have on the centrifugal-acceleration modes 
predicted by earlier analyses? 

(b) What "acoustic" modes are predicted due to fluid compressibility? 


The answers provided are as follows: 

(a) Compressibility has a negligible influence on centrifugal-acceleration modes. 
The result with and without compressibility are basically the same. 

(b) Fluid compressibility yields acoustic modes, with the lowest mode appearing at 
about twelve times running speed. This mode would be excited by rotor 
precession. 


Concerning pressure oscillations within the leakage annulus, Nhai used a Passing pre^ure 
wave at the pump impeller exit for excitation with n diametral modes. Following Bolleter 
(1988), n=n, (impeller vanes) - n 2 (diffuser vanes) = 24-13=11. Nhai s analysis predicts s^^arp 
peaks at precessional frequencies which are 6.5 and 7.8 times running spee • ese 
give predicted amplification factors from impeller discharge to the exit seal of the annulus 


(inlet seal of the impeller) of 7.5 and 17.8. Either mode could reasonably explain the internal 
melting observed in the HPFTP seals. The difficulty is that no excitation frequency is 
predicted by existing theories near these resonant frequencies. Bolleter predicts excitation 
frequencies at multiples of Q = n, ©/n where © is the running speed. Table 2 (page 38) of 
Nhai’s thesis shows no excitation frequencies at or near the resonant frequencies. 


To be brief, the present analysis only provides a plausible explanation for the melting-seal 
situation, if the impeller provides the required excitation frequencies, and testing would be 
necessary to confirm their presence or absence. 
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Introduction 

Figure 1 illustrates an impeller stage of a multistage cen- 
trifugal pump. Leakage along the front side of the impeller 
from impeller discharge to inlet, is restncted by a wear-nng 
seal, while leakage along the back side is restricted by ei 
an interstage seal or a balance-discharge seal. The axial thrust 
on the impeller is obviously of interest for structural mtegny 
of the pump and several investigators have presented analyses 
and test results for the thrust versus various impeller param- 
eters; Thomae and Stucki (1970), Lobanoff and Ross (1985). 
Impellers are sometimes used directly as thrust-balancing e 
ments; e.g., the main impeller of the Space Shuttle Mam Engine 
(SSME) High Pressure Oxygen Turbopump (HPOTP) t 
double-suction-entry design, and orifices at the inlet and exit 
of the leakage path create the principal axial-thrust-balance 
for the rotor. The leakage-path for the back side of the High 

Contributed by the Technical Committee on Vibration and Sound for put^ 
lication in the Journal of Vibration and Acoustics. Manuscript received 
January 1990. 


Pressure Fuel Turbopump (HPFTP) is also used for axial thrust 

^From an axial vibration viewpoint, the change in the axial 
thrust of an impeller which is used for axial-thrust balance is 
normally modeled by a stiffness and damping coefficient, i.e., 
F z =-KZ r -CZ r <» 

where Z is the axial change in position. The stiffness coefficient 
K ” fife tool slope in the thrust versus axial position curve 
and is nominally constant around the equilibrium position. In 
the SSME turbopumps, the stiffness of the pump housing 1 
used in seriS S the slope of the thrust-axial-posit.on curve 
to calculate *. Damping primarily arises due to flow throug 

"*« iSTaSmr.ion involved in .he 
(1) is that the natural frequencies of the fluid m the le * k ^ 
path are much higher than the pump’s running speedorthe 
axial vibration frequencies of the rotor. From a conventiona 





point R„Z, to an arbitrary point R,Z is denoted by S and 
defined by 

5 - hj" (§)’■'“ ' kJ 1+ (iV“ <Z) 

In the equations which follow, the path coordinate S and 
angular coordinate B are used as independent spatial variables. 
The coordinates Z,R defining the impeller surface are expressed 
as parametric functions of S, i.e., Z(S), #(S). The length o 
the leakage path along the impeller face is defined by 


[ Z ‘ +L ( dR\ 

Lj J z < 1 \dZ/ 


du 


( 3 ) 


Trigonometric functions of the angle 7, illustrated in Fig. 
3, are defined as follows 

dR dZ dR ... 

tanT= ■dZ• C0Sy= dS’ smy= ~ds () 

The clearance between the impeller and the housing is denoted 
as H(S,6,t), with the time dependency introduced by impeller 
motion. In the centered position, the clearance function de- 
pends only on S and is denoted by H 0 (S). Displacement of the 
impeller in the Z direction by the differential A yields 

dR 

= // 0 (S) - A(0sin7 = H 0 (S) + A(/>^ (5) 

Governing Equations 

Returning to Fig. 2, the path coordinate S and circumfer- 
ential coordinate RB are used to locate a fluid differential 
element of thickness H{S,6,t). From Childs (1987), the con- 
tinuity equation can be stated 

where U, and U e are the path and circumferential bulk-velocity 
components, respectively. Also from Childs (1987), the path 
and circumferential momentum equations are stated 

dP (4 dR 

- H Ts=- pH TTs +T ” +T " + 


ordinates (RB, S) and time, t. They are conveniently 
nondimensionalized by introducing the following variables 

u s = U s /V h u e = Ug/Rpi, P = P/p yl > 
h = H/C„ s = S/L s , r= R/R, (8) 

T = ut, b=V/R/*, T=L s /V, 

The objective of the present analysis is to examine the changes 
in (u„ u», p) due to changes in the clearance function h{B,s,t) 
caused by small axial motion of the impeller within its housing. 
To this end, the governing equations are expanded in the per- 
turbation variables 

u s =u 50 + eu sU h = h 0 + eh x 
u e = ue 0 + * u 6u P = Po + *P\ ^ 

where e = e/Cj is the perturbation parameter. The following 
equations result: 


<4 = ° 


Zeroth-Order Equations 

(a) Path-Momentum Equation 

dp 0 \(dr\ul„ | [ (°r+°A L — 

ds r\ds ) ft 2 L V 2 / K ds r ds 

(b) Circumferential-Momentum Equation 

+ [aXu<jo -r) + cju,o 1/2 = 0 

ds r ds 

(c) Continuity Equation 


(10a) 


(10ft) 


rh a u s0 = 1 


(10c) 


du,. 




HdP U (W* dU,^ JUe 

-!~d8 =Te ’ +T ° r Pl \dF SB R 3S Us R as) 

Following Hirs’ approach (1973), the wall shear-stress defi- 
nitions in these equations can be stated 

ms + 1 

r^pUlRTV+iUe/U,)^ 


T„=y/i^n 1 + We ~ RWU ,) 1 ) 


mr+ l 

1*1” (6) 


ms* 1 


T fc -y pU ) u 3 RT s U+WUf) 2 


T» r =jpU s (U,-Ro>)R7 r l 1 + m-Ru)/U:V) 


where 


R s = lHU s /v 


( 7 ) 


The continuity equation has been used to eliminate ^ from 

equation (10a). The momentum equations define the pressure 
and velocity distributions for a centered impeller position . They 
are coupled and nonlinear and must be solved iteratively. The 
initial condition for «„ o (0) is obtained from the exit flow con- 
dition of the impeller. The inlet and discharge pressure of the 
impeller are known and serve, respectively, as the exit (P t ) and 
supply ( ) pressures for the leakage flow along the impeller 
face. The inlet condition for p B is obtained from the inlet 
relationship 

P J -P o (0,M = Pd + *)C4(0,M/2 (1 1) 

From this relationship, the zeroth-order pressure relationship 
is 

p 0 ( 0) = P/pVf - (1 + 0«4( 0)/2 (12) 

The wear-ring seal at the leakage-path exit also provides a 
restriction, yielding a relationship of the form 

P(L s ,B,t)-P r =%C a ML s ,B,t) (13) 

First-Order Equations 

(a) Path-Momentum Equation 

^ + UeiA* + u sl A is +[uT^ + uT^^ + »so-£: “Mu 

(14fl) 

(ft) Circumferential-Momentum Equation 
ft^ 5 - + W»1 A-y) + U,\Am 

(c) Continuity Equation 

du s] uTdu„ (\ dr 1 dh 0 \ 


(14ft) 


Ai« JO dh 


h\ 


dh 0 1/ 3V i\ 

r 98 w7 aJ 


Nondimensionalization and Perturbation Analysis 

The governing equations define the bulk-flow velocity com- 
ponents ( U s , U e ) and the pressure P as a function of the co- 


(14c) 




Most of the parameters of these equations are defined in Ap- 
pendix A. The quantities a , and a r are defined by 

<j s =(L/H 0 )\„ a, = (£/ H 0 )\ r ( 15 ) 

where X, and X r are dimensionless stator and rotor friction 
factors defined by 


\ = nsR™[l+(u eo /bu so f] 2 

mr + | 

X, = nrR™ ( 1 + l(u eo - r)/bu so ] 2 ) ~ 

From equations (5) and (8), the perturbation clearance func- 
tion is 



Hence, 

dhi ^5 (Rj\(fr dh ] 
ds ~e\Ljds 2 ’ ~dd =0 

First-Order Equation Solutions. For axial excitation, the 6 
partial derivatives of equations (14) are eliminated. The time 
variation can be eliminated by assuming the following har- 
monic solutions for the clearance excitation 

& = & 0 e i °' = 5 0 e i/r ( 17 ) 

and dependent variables 

Pi=Pie J/r , u sl = u sl e ifr , ut l ^u n e if ' ( 18 ) 

where 

f=Q/w, <j)t = T (19) 

Note that p x , u s[ , and u 91 are now complex variables. Sub- 
stitution into equations (14) yields 


where 



IA] = 


B 0 

Ayg/u s0 (A 2e +JQT)/u s , 

s-UsoB+jQT At, 


SIT dr 




A dr . 0 Tdr 

Au di* u " F '* ,u “T.* 


and 


■^!ho(fL_±^odr\ 1 dr I _dh 0 

h 0 yds 2 h 0 ds ds) * r ds + h Q ds 


( 20 ) 


( 21 ) 


( 22 ) 


( 23 ) 


The following three boundary conditions are specified for the 
solution of equation (20): 

(a) The entrance-perturbation, circumferential velocity is 
zero, i.e.. 


«»i(0) = 0 (24 a) 

(b) The entrance loss at the seal entrance is defined by equa- 
tion (II), and the corresponding perturbation-variable rela- 
tionship is 

Pt(0)= -(1 + $)u si (0) (24b) 

(c) The relationship at the exit is provided by equation (13) 
and yields the following perturbation relationship 
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Table 1 Zeroth -order-solution results; C,= 3.5 mm 


«#o(0) 

c* 

m {kgf$te} 
F„ x 10“ 4 (AT) 


0.5 

1.655 

4.784 

2.25 


0.6 

1.600 

4.301 

2.15 


0.7 

1.731 

3.073 

2.04 


Table 2 Asymptotic coefficients from Figs. 5 and 6 


UtfO 

K. x to -4 

M 0 

C„ x 10 3 


N/m 

Kt 

Kg 

0.5 

.648 

5.61 

.809 

0.6 

6.06 

5.63 

.911 

0.7 

13.6 

5.68 

1.05 


P.0) = CW 1)MD (24c) 

The value for C de depends on the wear-ring seal geometry. 
Solution of equation (20) for the boundary conditions of equa- 
tions (24) is relatively straightforward using a transition-matrix 
approach (Meirovitch. 1986). The solution can be stated 

^5 A C/T) (25) 


/s \ ! 

Tic+XfiT 

>,-(=) 

Ar+i/if ! 




- ? M 


(A) 


Reaction Force. From Fig. 4, the axial differential force 
component on a differential-impeller surface area can be stated 

dF z = - (Psiny + r n cosy)RdddS (26) 

The zeroth and first-order force components are obtained by 
substituting for P and t„ in equation (26). From equation (6), 
the perturbation shear stresses can be stated 

V, = P y?(B sl u sl + B a u e{ + B ti hJ (27) 
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The coefficients of these equations are defined in the appendix. 
The zeroth and first-order forces are defined by 


f 

r ,,4r /, 

l\dz 

Jo 

ns** 0 

l) dr Trs0 


rds 


and 

where 


F zi ~ Fz\k +jF Z \ c = F 0 (f k +jff c )5 0 
F 0 =irR}AP=TRfC d pV} 


4k'[4 


(28) 

(29) 

(30) 


rds (31) 
rds 



Numerical Results. Figure 4 illustrates the pump-impeller 
and shroud geometry used by Bolleter et ai. (1987) in their test 
program for radial force coefficients. Their pump uses a vaned 
diffuser. Their tests were at best efficiency point (BEP) with 
the pump running at 2000 rpm, while developing 68m of head 
and 130 I/sec of flow rate. The impeller has seven blades and 
an ^impeller exit angle of 22.5 deg. The test fluid is water at 
80 F. For the present study, AP across the impeller is assumed 
to be 70 percent of the total head rise of the stage. Based on 
pitot-tube measurements, impeller exit tangential velocity is 
about 50 percent of the impeller discharge surface velocity; 
hence, w*,(0) = 0.5. 

Both walls of the annulus are assumed to be smooth and 
represented by Yamada’s (1962) test data; mr = ms= -0.25, 
nr = ns= 0.079. The inlet loss for the annulus, £, is assumed 
to be 0.1. The discharge coefficient for the seal is calculated 
iteratively as follows. With an assumed C de , equations (10), 
(11), and (12) were used to calculate the leakage through the 
impeller annulus and the pressure and tangential-velocity up- 
stream of the seal. The seal is then analyzed (with the same 
equations) using the calculated seal inlet pressure and tangen- 
tial velocity to determine leakage and C*. The iteration con- 
tinues until the leakage predictions for the exit seal and the 
impeller annulus agree. Table 1 provides zeroth-order solu- 
tions. 

Figures 5 and 6 illustrate/* and/, versus /for 1^(0) = 0.5, 
0.6, and 0.7. The «#o(0) = 0.5 curves are comparatively smooth; 
the w*o(0) = 0.6, 0.7 curves show evidence of fluid resonances 
similar to those obtained earlier by Childs (1987). Specifically, 
in the absence of fluid resonances, the expected results for/* 
would be a parabola without the predicted fluctuation in the 
neighborhood of /= 0.25. Further, the expected result for / 
would be a constant without the low-frequency fluctuations. 

The functions /*, / are nondimensionalized frequency-de- 
pendent stiffness and damping coefficients. To develop a phys- 
ical model for the axial reaction forces defined by these curves, 
the /*(/) curves will be reviewed first. All of the /* curves 
demonstrate a quadratic asymptote with the following fre- 
quency-domain model 

fko<J) = -K+Mf, (32) 

which implies the time-domain reaction-force model 

F* 1J =-(/ir ff Z, + M a Z r ) (33) 

where 

K q = KF 0 /C h M a = MF 0 /C?j? (34) 

The physical coefficients obtained from a curve fit of the 
asymptotic solutions in Fig. 5 yield the physical coefficients 
of Table 2. 



The stiffness values are negligibly small in comparison to 
the axial stiffness which would center the pump rotor either 
through a thrust bearing or balance-piston arrangements. The 
mass coefficient is small but appreciable in comparison to the 
impeller mass. This “added mass” contribution to the impeller 
rotor has not previously been accounted for in axial vibration 
analyses of pumps. 

The / curves of Fig. 6 can be viewed as nondimensionalized 
frequency-dependent damping coefficients. The asymptotic 
behavior of these curves (for high values of f) shows an un- 
expected, approximately-linear, increase in damping with in- 
creasing frequency. The asymptotic results are curvefitted by 
the linear model 

fc<f) = C 0 -Cf (35) 

The linear dependence of / on C yields the time-domain re- 
action-force model 

F ca = - CtfiaZr* = I Z/Z, I 1/2 (36) 

where 

C^CFJC^ 1 

Values of C a are provided in Table 2. Physically meaningful 
values for damping are obtained by multiplying C a by Q, e.g., 
the nominal damping value for Q = w = 209 rd/sec is about 209 
N sec/ m (1.2 lb sec/in) which is significant. 





The physical example used to demonstrate the solution of 
equation (39) is a double-suction single-stage pump with the 
impeller of Fig. 4 and the following physical data: 

M zo = 100 Arg, f = 0.05, u = 2000 rpm = 209 rdf sec 

With a double-suction design, the reaction force from the 
impeller face is doubled. Figures 7 and 8 illustrate solutions 
to equation (39) for ** = 0.5 and 1.5 times the running speed 
for no impeller forces and impeller forces corresponding to 
Ww (0) = 0.5, 0.6, and 0.7. The magnitude and phase of * are 
presented and demonstrate that the impeller forces provide a 
modest amount of additional damping and slightly displace 
the peak-amplitude-frequency location. Although not pre- 
sented here, the impeller forces did cause significant changes 
in the steady state response when *„ was reduced to unreal- 
istically low values on the order of 0.1. 

Summary and Conclusions 

An analysis has been developed and predictions presented 
for the axial forces developed on a pump impeller shroud. The 
force coefficients are reduced to frequency-dependent stiffness 
and damping coefficients. Fluid resonances are evident in these 
results comparable to the rotordynamic-coefficient analysis 
results of Childs (1987). The asymptotic expansion of these 
coefficients yield negligible stiffness values and small but ap- 
preciable added mass and damping coefficients. 

Frequency-response analysis of a double-suction single-stage 
pump impeller indicates that the impeller shroud forces provide 
a little additional damping and can modestly move the damped 
natural frequency of the system for pumps with axial natural 
frequencies on the order of 0.5 to 1 .5 times the running speed. 
The impeller-shroud .forces can only become really significant 
for natural frequencies that are much lower than the running 
speed, viz., one tenth of the running speed or lower. 
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Modeling Annular Forces for Axial-Vibration Analysis. The 
transient model for axial vibrations, including the impeller 
forces, can be stated 

+ (C z0 - F c )Z r + (K„ - F k )Z r = F'(» (37) 

where M. a , C zo , and K zo are the nominal mass, damping, and 
stiffness coefficients, F e (t) is the external excitation force, and 

Ftf a ) =f'tf a )F 0 /C l w, F k (f„) =/*0 QFJC t (38) 

/ fl =n fl /w 

The model of equation (2) has moved the impeller-shroud 
reaction force from the right-hand side of the equation to >the 
left-hand side of the equation. The force is now modeled as 
frequency-dependent stiffness and damping coefficients, which 
combine directly with the nominal, mechanical, stiffness and 
damping coefficients. 

For an external harmonic excitation force of the form 
the steady-state solution = is defined by 
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appendix a 

Perturbation Coefficients 

A i, = [of 1 - ms) + of 1 - mr)] u 2 so /2h 0 


a — \a\^ -A /(F*/ K zo ) 

= a ) \/ ( [*£( 1 - /*) -fi] +i/w /J (2f - / c ij„) ] 


where 


u 2 n = K zo /M zo \ 2 (o>„ = C zo /M zo 


(39) 


(40) 

, du so 

Ais ~ ds 


r ds 


- [(1 + mr)afl 0 (u Bo - r) + (1 + ms)ofi x u Bo ]/2 
2Aj 0 = « J0 [(1 - mr)(u 6o - r)o r + (1 - ms)u Bo o s \/h 0 
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2 A u = u so {a r + a s ) + a^mr + 1 ){u 6o - r)0 g 

+ afms+l)ue o 0 > +2^^f 
r as 

2/4 = aXu« 0 -r)[mr-( 1 + mr)0 g (u eo - r)/u so ] 

+ oMclms-V + ms)0 t u $o /u so ] 
00 = («fe - r )/b 1 u w { 1 + [(u 0o - r)/buj 2 ) 

= W( 0 /fe J w J0 [l + {Ug Q / bu s0 ) 2 } 


th> Perturbation Coefficients 

B ey =X^1 + mr)(u eo -r)[l - 0 „(u eo - r)/u s0 ]/2b 
B n = Kl“s 0 + (1 + mr)(u 6o - r)0 o ]/2b 
Be 3 = X/ nr(Ue Q - r)uJ2bh 0 
t„ Perturbation Coefficients 

B s i = X r [( 2 + mr)u so - ( 1 + mr)/3 0 (u„ 0 - r)/2] 
= X r (l + mr)P 0 u so /2 
0 s i = \snn?J2h 0 




^ o 

^ T 
\D H- 
vD I — 1 
(- 1 P- 
CT W 




| 



I I 


N 


i 



! 




Reprinted from April 1991, Vol. 113, Journal of Vibration and Acoustics 




D. W. Childs 

Turbomachinery Laboratories, 
Mechanical Engineering Department, 
Texas A&M University, 
College Station, TX 77843 


Centrifugal-Acceleration Modes for 
Incompressible Fluid in the 
Leakage Annulus Between a 
Shrouded Pump Impeller and Its 
Housing 

An analysis is presented for the perturbed flow in the leakage path between a 
shrouded-pump impeller and its housing. A bulk-flow model is used for the analysis 
consisting of the path-momentum, circumferential-momentum, and continuity equa- 
tions. Shear stress at the impeller and housing surfaces are modeled according to 
Hirs’ turbulent lubrication model. The governing equations have been used earlier 
to examine rotordynamic reaction forces developed by lateral and axial impeller 
motion. 

A perturbation expansion of the governing equations in the eccentricity ratio 
yields a set of zeroth and first-order governing equations , The zeroth-order equations 
define the leakage rate, and the velocity and pressure distributions for a centered 
impeller position. The first-order equations define the perturbations in the velocity 
and pressure distributions due to axial or lateral motion of the impeller . Prior analyses 
by the author of the perturbation equation have examined the reaction forces on 
the shroud due to rotor motion . These analyses have produced “resonance” phe- 
nomena associated with the centrifugal-acceleration body forces in the fluid field. 

In the present analysis, an algorithm is developed and demonstrated for calculating 
the complex eigenvalues and eigenvectors associated with these resonances. First - 
and second-natural-frequency eigensolutions are presented for mode shapes corre- 
sponding to lateral excitation. First-natural-frequency eigensolutions are also pre- 
sented for mode shapes corresponding to axial excitation . 


Introduction 

Figure 1 illustrates an impeller stage of a multistage cen- 
trifugal pump. Leakage along the front side of the impeller, 
from impeller discharge to inlet, is restricted by a wear-ring 
seal, while leakage along the back side is restricted by either 
an interstage seal or a balance-piston discharge seal. The pres- 
ent analysis considers perturbed flow in the leakage paths be- 
tween the impeller shroud surface and its housing. 

Prior analyses by the author of those annulli have been 
concerned with lateral (1987, 1989) and axial (1990) reaction 
forces developed by the impeller shrouds as a consequence of 
impeller motion. These analyses have been based on “bulk- 
flow” models which neglect the variation in the dependent 
variables across the fluid film. The model consists of the path 
and circumferential momentum equations and the continuity 
equations. 

The analyses cited have yielded force and moment coeffi- 
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cients due to impeller motion but have also predicted “reso- 
nance” phenomena, which are caused by the centrifugal- 
acceleration body forces present in the path momentum equa- 
tions. In the present analysis, an algorithm is developed and 
demonstrated for calculating the complex eigenvalues and ei- 
genvectors associated with the fluid resonances. 

Geometry and Kinematics 

Figure 1 illustrates the annular leakage paths along the front 
and back sides of a typical shrouded impeller of a multistage 
centrifugal pump. The present discussion concentrates on the 
flow and pressure fields within the forward annulus; however, 
the analysis also applies to the rear annulus. As illustrated in 
Fig. 2, the outer surface of the impeller is a surface of revo- 
lution formed by rotating the curve R = R(Z) about the Z axis. 

A point on the surface may be located by the coordinates Z, 

R(Z),e. The length along the curve R(Z) from the initial point 
R h Zj to an arbitrary point R,Z is denoted by S and defined 
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In the equations which follow, the path coordinate S and 
angular coordinate 0 are used as independent spatial variables. 



The coordinates Z t R defining the impeller surface are expressed 
as parametric functions of S, i.e., Z(S),/?(S). The length of 
the leakage path along the impeller face is defined by 



( 2 ) 


Trigonometric functions of the angle 7, illustrated in figure 2, 
are defined as follows 


tany= - 


dR 

dZ' 


dZ 

COSy = — , 


sin7 = - 


dR 

dS 


(3) 


The clearance between the impeller and the housing is de- 
noted as with the time dependency introduced by 

impeller motion. In the centered position, the clearance func- 
tion depends only on S and is denoted by H 0 (S). Displacement 


Nomenclature 


AietAtftAtf 

Q, 

c, 

C r 

f- 0/u) 


fkJc 


frq,fdq 


= coefficients intro- 
duced in equation 
(14) and defined in 
the appendix 
= discharge coefficient 
for the exit wear- 
ring seal, introduced 
in equation (13) 

= initial (s = 0) radial 
clearance (L) 

= exit seal clearance 

(L) 

= dimensionless exci- 
tation frequency 
= nondimensional 
stiffness and damp- 
ing coefficients for 
the impeller corre- 
sponding to axial 
motion at the non- 
dimensional fre- 
quency / = 0 /oj 
= nondimensional, ra- 
dial, and circumfer- 
ential impeller- force 
coefficients corre- 
sponding to a circu- 
lar orbit at the 


nondimensional fre- 
quency^ o 7« 

h-H/Q = nondimensionalized 
clearance 

H — clearance between 
impeller shroud and 
housing (L) ^ :r . 


L s = leakage-path length, 
defined by equation 
(2), (L) 

p-P/pV) = nondimensionalized 
static fluid pressure 

P = fluid static pressure 
(F/L 2 ) 

r=R/Ri ~ nondimensionalized 
radial coordinate 

R = radial coordinate 

(L) 

Ri = initial (y = 6) radius 

( L ) 

R s = 1HU s /p - path-velocity Rey- 
nolds number 

s = S/L s - nondimensionalized 
path length 

5 = path coordinate in- 
troduced in equa- 
tion (1), (L) 

T=L S /V i = representative 

transit time for fluid 
traversing the leak- 
age path (7) 

u s =U s /Vj - nondimensionalized 
path fluid velocity 
u 9 = Ue/Rfu = nondimensionalized 
circumferential fluid 
velocity 

V i = initial (5 = 0) rotor 
velocity 

a = a+yw rf = dimensionless com- 
plex eigenvalue for 
fluid mode 


e = perturbation coeffi- 
cient 

w = pump running speed 
(T- 1 ) 

Q = excitation frequency 
(T~ ] ) 

p - fluid density (A i/V) 

6 = circumferential co- 
ordinate 

o„o 5 = normalized friction 
factors, defined by 
equation (15) 

7 — (jit * nondimensionalized 
time 

£ = entrance-loss coeffi- 
cient introduced in 
equation (11) 

A - axial impeller dis- 
placement ( L ) 

<r=fa n - dimensionless real 
part of complex ei- 
genvalue for fluid 
mode 

{* = damping factor 

from complex eigen- 
value 

= dimensionless imagi- 
nary part of com- 
plex eigenvalue for 
fluid mode (damped 
natural frequency) 
u n = dimensionless un- 
damped natural fre- 
quency from 
complex eigenvalue 
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of the impeller in the X and Y directions obviously causes a 
change in the clearance function, as does a change in the axial 
position defined by A(/). For small displacements and rotations 
of the impeller, the clearance function can be stated 

= H 0 (s) - Xcosd - YsinO - Asiny (4) 

Observe in this equation that H 0 and siny are solely functions 
of 5, while X , Y t and A are functions only of /. 

The clearance function provides the excitation for the re- 
action forces developed in earlier analyses. Its importance for 
the present study concerns the nature of the assumed solution 
in seeking eigen solutions. 


Governing Equations 

Returning to Fig. 2, the path coordinate S and circumfer- 
ential coordinate R8 are used to locate a fluid differential 
element of thickness H(S,6 t t). From Childs (1987), the con- 
tinuity equation can be stated 




(5a) 


where U s and U e are the path and circumferential bulk-velocity 
components, respectively. The path and circumferential mo- 
mentum equations are stated 


-a£-r. 

3S 


„lA dR „ 

+ T ”- pH -£li + pH 


dUs 3C/, Ug dU ,.. \ 
,3t 36 R dS 7 


(5b) 


H3P 
R 36 


Tds + T 6r + pH 


/3Ue dUeUe 
V 3 1 + de r 


dUe U e U s dR\ 
dS 1 R 35/ 


(5c) 


Following Hirs* approach (1973), the wall shear-stress defi- 
nitions in these equations can be stated 


ms ♦ 1 

r a =^pUlRril +{U t /U t fP~ 


Tsr=j(>t£R7 r {l + [(U,-Ru)/U s ] 2 ) 2 (6) 


ms + 1 

r* = tfpUeU'RT’V + (Ug/U,) 2 ~ 


u s = u so + eu sU h = h 0 + eh t 

u e = u eo + eu t ,, p=p 0 + tPi (9) 

were e = e/C, is the perturbation parameter. The following 
equations result: 


Zeroth -Order Equations 

(a) Path-Momentum Equations 


(dr\i4 0 T 

/o r + a A 1 dh 0 

1 dr 

\ds) b 2 L 

\ 2 ) h 0 ds 

r ds_ 


(b) Circumferential-Momentum Equation 

— + — % + - r) + ffjW« 0 ]/2 = 0 

ds r ds 

(c) Continuity Equation 

rh 0 u S0 = 1 


(10a) 


(106) 

(10c) 


The quantities a, and a r are defined by 

a, = (L/H 0 )\ s , a, = {L s /H 0 )\ (11) 

where \ and X, are dimensionless stator and rotor friction 
factors defined by 


Xj = + (ugo/buj 2 ] 

\ r =nrRZ r [l + [(u<, 0 -r)/bu S0 ] 7 \ 


mr+ 1 


du„ 


The continuity equation has been used to eliminate from 

equation (10a). The momentum equations define the pressure 
and velocity distributions for a centered impeller position. They 
are coupled and nonlinear and must be solved iteratively. The 
initial condition for w flo (0) is obtained from the exit flow con- 
dition of the impeller. The inlet and discharge pressure of the 
impeller are known and serve, respectively, as the exit ( P e ) and 
supply ( P s ) pressures for the leakage flow along the impeller 
face. The inlet condition for p 0 is obtained from the inlet 
relationship 

P, - P o (0,e,t) = p( 1 + t)U] o (0,8,t)/2 


From this relationship, the zeroth-order pressure relationship 
is 


p o (0) = P/p^-(l + i)id(0)/2 (12) 


mr+ 1 

r*= jpUAU,-Ru>)R7 r { 1 + l(U'-Ro>)/U s ] 2 \ 1 


The wear-ring seal at the leakage-path exit also provides a 
restriction, yielding a relationship of the form 


where 

R s = ZHU/v 


(7) 


PiLJj) -P t = ^Q,C?(T„0, 0 


(13) 


Nondimensionalization and Perturbation Analysis 

The governing equations define the bulk- flow velocity com- 
ponents (£/„ U e ) and the pressure P as a function of the co- 
ordinates (R6 f S) and time, t . They are conveniently 
nondimensionalized by introducing the following variables 

t h=U s /V i , u e =U,/R l u, P = P/py} 
h = H/C), s=S/l s , r=R/Rj (8) 

r = ut, b=V/R,u>, T=LJV, 

The objective of the present analysis is to examine the changes 
in (w„ u 9t p) due to changes in the clearance function h(6 f s , 
/) caused by small motion of the impeller within its housing. 
To this end, the governing equations are expanded in the per- 
turbation variables 



f * w m i 


Annn ^nn- 1 -110 1 on 



First-Order Equations 

(a) Path-Momentum Equation 
dp\ 

— + U ei A 2 s + U s \Ais + 


. J*eo du s\ . Mi 

0)T-— + wT — +Mj„ , 

L St r dd os . 


- 

(14a) 


(6) Circumferential-Momentum Equation 


,L S 1 3pi . . 

^ + + u j1^3(9 


wr— +CJP ^“ + w J0 — = Mi<? (146) 

ot r oo os 


(c) Continuity Equation 
du s i 
ds 


<j)T du e t 1 t/AA 

r dB + Usl \r ds + h 0 ds ) 

h x u s0 dh 0 1/ 3A, u 6t 

^~n~ 9 ~r-~r( w JO — -f o>7^- 
/to ds 6 0 \ ds r 


a6, 

aa 


+ o)T- 


dr) 


(14c) 


New coefficients in these equations are defined in Appendix 
A. 


Eigen Solutions Corresponding to Lateral Excitation 

The first-order equations (14) define the first-order pertur- 
bations u 5 [ (s t $ t T) t «(n(5,a t r), and p(s y d,T) resulting from the 
perturbed clearance function 6|(r). From equations (4) and 
(8), h x can be stated 

eh i = 6ic($,t)cos0 + h ls (s,T)sind (15) 

The 0 dependency of the dependent variables is eliminated by 
assuming the comparable, separation-of-variable, solution for- 
mat 

u s \ = u s]c cosd + u sls s\nB 

u ei = u^cCosO + u 9ls smd (16) 

Pi - + PijSinfl 

Substituting into equations (14) and equating like coefficients 
of cos# and sinfl yields six equations in the independent vari- 
ables 5,r. By introducing the complex variables 

!4ri = «,ic +Mi„ Hex = ueic +M\ S 

E\=Pic+jPis. h\ =h ic +jh ls , (17) 

these six real equations are reduced to the following three 
complex equations in s and r. 

^ + »i> 4* + «,,/!*+ -juT^u^+uJ^ =h,A u 

(18a) 


-Mu (186) 


-A 0jr) £i + UnAu + UjnA } , 

+ [ w7 l? - JuT T-' ,+ u “^ L _ 

du„ uT (\ dr 1 dh\ 

lr-^*' + ±'Krds + h 0 l?) 

\}hodho ,vTu»o\ 

-‘ [ ds ^67 r J h 0 ds h a br C 


The time dependency of equations (18) is eliminated by as- 
suming a harmonic solution of the form 

6, = A l0 e^, a,, = u,, = £, =p,e^ T ; /= Q/w (19) 






r ■ ‘h* 

Fig. A Nondlmenslonal rotordynamfc force coefficients: (e) tangential- 
force coefficient, (0) radial-force coefficients 


The assumed clearance function h x (t) correspond to circular 
orbital motion of the impeller at the precession frequency Q. 
Substitution from equation (19) into equation (18) yields the 
three complex ordinary differential equations in s of the form 



The parameters q 0t g u g 2 » and S 3 on the right hand side are 
defined by Childs (1987). The right hand side is important for 
the forced solution but irrelevant for eigenanalysis, which con- 
cerns the homogeneous solution. 

The following three boundary conditions are specified for 
the solution of equation (20): 

(a) The entrance-perturbation, circumferential velocity is 
zero, i.e., 

uex( 0) = 0 (21a) 

(b) The entrance loss at the seal entrance is defined by equa- 
tion (12), and the corresponding perturbation- variable rela- 
tionship is 

Pt(0)=-(l+£)«,i(0) (21 b) 

(c) The relationship at the exit is provided by equation (13) 
and yields the following perturbation relationship 

p l (l) = C d ^ J0 (l)u Jl (l) (21c) 

The value for C& depends on the wear-ring seal geometry. 

Figure 4 illustrates the radial and tangential reaction-force 
components for solutions to equation (20) from Childs (1989). 
The three solutions for each frame of the figure correspond 
to three inlet tangential velocity ratios w w (0). If the centrifugal 
acceleration term is dropped from the analysis, the f 9q curve 
becomes a straight line and the curve becomes a parabola, 
which are the expected solutions for liquid, constant-radius, 
annular seals; Childs and Kim (1985). The sharp deviation 
from the expected linear and parabolic solution forms for these 
curves in the “ resonance*’ phenomenon cited in the introduc- 
tion, and it is eliminated if the centrifugal acceleration term 
is dropped. The question of interest here is: “If the solutions 



of Fig. 4 arise from a ‘resonance* of the fluid, how are the 
associated eigenvalues and eigenvectors to be calculated, and 
what do the mode shapes look like?” 

The first part of the question is answered by reviewing the 
solution approach to equation (20) which was used in gener- 
ating the results of Fig 4. Following normal transition-matrix 
approaches (Meirovitch, 1986), the solution to equation (20) 
can be stated 

(•*(*)) = [*(s)J ( Wo ) + { V(s ) ) ; (*v) r = (Mji, u„„ p { ) (22) 

where [$] is the transition matrix, (w 0 ) is the vector of initial 
conditions, and (u) is the particular solution obtained for zero 
initial conditions. [$] is obtained by solving the homogeneous 
version of equation (20) three times for the initial conditions 
(1,0,0), (0,1,0), and (0,0,1). One of the three initial conditions 
of equation (22) is given; viz., w<n(0) = 0. Equation (216) pro- 
vides one relationship between the remaining unknown initial 
conditions pj(0) and u 5l (0). The final relationship between 
these variables is obtained by evaluating equation (22) at 5 = 1 
to obtain 

isi(iy = *nO)«si( 0 ) + *1 3 (1)P|(0) + »,( 1 ) 

Pi( 1 ) = * 3 |< 1 )»,|( 0 ) + *»(l)p,(Q) + u 3 (l), ( 23 ) 

and substituting into equation (21c). The resulting set of equa- 
tions for the unknown initial conditions can be stated 

(i + O i If M0)] 

i) ^(d-c^dJUko)) 

0 

-^(O + C^UjotOvid) 

For forced response, the right-hand side is nonzero, and the 
missing initial conditions are obtained by inverting the coef- 
ficient matrix. 

The resonance peaks in Fig. 4 appear when / is near the 
imaginary part of the complex eigenvalue. To obtain the ei- 
genvalue, the harmonic solution of equation (19) is replaced 
by the general solution format 

Kfi = «n« aT - H9\ = u n er t £j =p,e aT 
where a is both a complex number and the desired eigenvalue. 
Substitution into equations (18) yields the following definition 
for the coefficient matrix of equation (20) 





Fig. 5 Nondlmenslonal axial fore# coal flclenta: (a) stiffnaaa force coaf- 
ficlant, (6) damping-force coefficient 




B 

Ais~ Wjofl+yTT 


_ &T 
* r 

(A-u+jTTl/ua 

A UT 

Ais+J WjO 




(25) 


where 

'-(?£*£*)• <“> 

Now, the differential equations and their homogeneous so- 
lutions are a function of the desired eigenvalue a, and the 
eigenvalue is found when the determinant of the coefficient 
matrix of equation (24) is zero. 

The following approach was used to solve for the complex 
eigenvalues: 

{a) An initial value, ao, is guessed based on the results of Fig. 
4; viz., ct o = 0+jf o where/, yields a peak on the f 9q curve. 

(b) The homogeneous version of equation (20) is solved to 
obtain [$], and the determinant 

De = ( 1 + ?)[* 33 <1) - CWl)] - * 31 <I) + C ae uAD (27) 
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Fig. 7(a) t/„(t) f Irat-natyraMrequency eigenvector corresponding to lat- 
eral excitation; C,*s 0.36mm 




Fig. 7(6) Ufi{s) flrst-naturaldrequency eigenvector corresponding to lat- 
eral excitation; C,= 0.36mm 


is calculated at a 0 and c*i = a 0 [l +5(1 +j) 1 to obtain De Q and 
De j; 5 is a small parameter. 

(c) The secant rule is then used to calculate the next estimate 
a 2 and Dei and subsequent a/s and De ' s until convergence is 
achieved. 

The eigenvalues are obtained by setting 

w,t(0)=l, p,(0) = (l+£), wai(0) = 0 

and evaluating equation (22) 

( w(s) }=[*(*,«) 1M0)} 

over «[0,1]. 


Eigen Solutions Corresponding to Axial Excitation 

For axial excitation, the perturbed clearance function h x {j) 
is defined from equations (4) and (8) to be 

A(T)/R\dr 




This excitation is not a function of 5, and neither are the 
associated forced-response solutions. Figure 5 illustrates real 
and imaginary reaction-force coefficients resulting from axial 
impeller excitation (Childs 1990). Note that “resonance’ phe- 
nomena are also present in these results. 

To obtain eigenvalues and eigenvectors, the following cor- 
responding assumed solution is substituted into equation (14) 

u sl = u sl e ar 1 Pi=Pi^ r » 


and the 6 derivatives are dropped 
equations 

. f 

a \ - / r . , 
T\u ei \ +[AJ 
ds / - \ 


to obtain the 



= 0 


homogeneous 




Fig. 7(c) p,(s) firsl-naturel-frequency eigenvector corresponding to lat- 
eral excitation; C r « 0.36mm 


214/ Vol. 113, APRIL 1991 


Transactions of the ASME 


Table 1 Zeroth-order-solutlon results; C,» 3.5mm, C r * 0.36mm 


tt#o(0) 

0.5 

0.6 

0.7 

C* 

1.655 

1.690 

1.731 

m (kg/ see) 

4.784 

4.391 

3.973 


Table 2 Lowest-natural-frequency eigenvalues corresponding to lat- 
eral excitation 


«#o(0) 

C r (mm) 

a 

“d 



0.5 

.36 

-.621 

.196 

.953 

.651 

0.6 

.36 

-.288 

.300 

.693 

.416 

0.7 

.36 

-.174 

.349 

.446 

.391 

0.5 

.72 

-.670 

.175 

.936 

.716 

0.6 

.72 

-.315 

.291 

.734 

.429 

.07 

.72 

-.188 

.345 

.480 

.392 



Table 3 Second-natural-frequency eigenvalues corresponding to lat- 
eral excitation; C,= 0.36mm 


u*o(0) 

u 

Wd 

f 

Un \ 

0.5 

-.787 

1.2 

.549 

1.43 

0.6 

-.344 

1.21 

.273 

1.26 

0.7 

-.214 

1.23 

.172 

1.25 




Fig. 8(a) </„(*) second-natural-frequency eigenvector corresponding to 
lateral excitation; C , « 0.36mm 


where 



r b 

0 

o' 

KJ- A v /u s o 

(Aw + awT^/Uso 

0 


A* 

0 


(28) 


Equation (21) continues to define the boundary conditions, 



Fig. 8(6) i> #1 (s) second-natural-frequency eigenvector corresponding to 
lateral excitation; C f « 0.36mm 




Fig. 8(c) p,(a) second-natural-frequency eigenvector corresponding to 
lateral excitation; C r - 0.38mm 
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Table 4 Firsl-natural-frequency eigenvalues corresponding to axial ex- 
citation 


Utfo (0) 

C r (mm) 

a 

Wd 

f 

Wn 

0.5 

.36 

-1.06 

.478 

.911 

1.16 

0.6 

.36 

-.584 

.396 

.878 

.706 

0.7 

.36 

-.400 

.364 

.740 

.541 

0.5 

.72 

-1.12 

.505 

.912 

1.23 

0.6 

.72 

-.619 

.418 

.829 

.746 

.07 

.72 

-.680 

.382 

.740 

.568 


and an eigenvalue is obtained when the determinant of equation 
(27) is zero. 

Numerical Results 

Introduction. Figure 6 illustrates the pump-impeller and 
shroud geometry used by Bolleter et al. (1987) in their test 
program for radial force coefficients. Their tests were at best 
efficiency point (BEP) with the pump running at 2000rpm, 
while developing 68m of head and 130//sec of flow rate. The 
impeller has seven blades and an impeller exit angle of 22.5 
deg. The test fluid is water at 80°F. For the present study, A P 
across the impeller is assumed to be 70 percent of the total 
head rise of the stage. Based on pitot-tube measurements, 
impeller exit tangential velocity is about 50 percent of the 
impeller discharge surface velocity; hence, w w = 0.5. 

Both walls of the annulus are assumed to be smooth and 
represented by Yamada's (1962) turbulence-coefficient test 
data: tnr = ms — -0,25, nr = ns = 0.079. The inlet loss for the 
annulus, £, is assumed to be 0.1. The discharge coefficient for 
the seal is calculated iteratively as follows. With an assumed 
C de , equations (10) through (13) are used to calculate the leak- 
age through the impeller annulus and the pressure and tan- 
gential-velocity upstream of the seal. The seal is then analyzed 
(with the same equations) using the calculated seal inlet pres- 
sure and tangential velocity to determine leakage and C de . The 
iteration continues until the leakage predictions for the exit 
seal and the impeller annulus agree. Table 1 below provides 
zeroth-order solutions. 

Eigen Solution Results Corresponding to Lateral Excitation. 

Table 2 cont ains the complex eigenvalues a = o + ju d - 
- + ju„yj\ - p for the lowest-frequency eigenvalue. Re- 

sults are provided for both nominal (C r = 0.36mm) and double 
radial seal clearances. Observe that the damped natural fre- 
quency <j) d ranges from about 20 percent to 30 percent of run- 
ning speed. Starting guesses for the eigenvalues were 
a 0 = 0+0.4 /, because the lowest-frequency peaks in Fig. 4(a) 
are around/=0.4. Observe that the calculated eigenvalues are 
entirely consistent with the results of Fig, 4(a); specifically, 
the solutions are stable and highly damped, the damping factor 
f decreases as ««(0) increases, and both and the / value 
corresponding to peaks in the f 9q curves decrease as u w (0) 
increases. The results are observed to be relatively insensitive 
to changes in the exit-seal clearance. 

Figure 7 contains the complex eigenvectors for the 
C r = 0.36mm. The u sl (s) vector has a peak amplitude at the 
inlet and exit with a small phase shift at the inlet. The w^i(s) 
magnitude and phase tend to increase steadily with increasing 
s. The Pi(s) magnitude increases steadily with increasing s until 
about s = 0.9 and then drops abruptly as it approaches the exit 
seal. The sharp change in u sX and p x between s = 0.9 and s= 1.0 
is caused by the sharply convergent flow path, illustrated in 
Fig. 5, as the flow approaches the exit seal. Note that the phase 
of a,! is near zero until it approaches the exit and then drops 
rapidly. The rapid phase shift at the exit is indicative of the 
large energy dissipation associated with the exit discharge seal. 

Table 3 contains eigenvalues for the second-natural-fre- 
quency eigenvector. Again, the damping factor and undamped 


natural frequency drop as « w ( 0) increases. Figure 8 illustrates 
the complex eigenvectors for these eigenvalues. Comparing 
u s] (s) in Figs. 7(a) and 8(a) shows similar mode shapes; how- 
ever, the peak amplitudes increase with increasing u m ( 0) for 
the first eigenvector but decrease with increasing w w (0) for the 
second eigenvector. The phase of Wotfs) increases for the first 
eigenvector with increasing s but decreases for the second ei- 
genvector. The Pi(s) amplitudes are similar for both the first 
and second eigenvectors; however, their phase behavior is quite 
different. 

Diligent searches revealed no eigenvalues with natural fre- 
quencies below the first eigenvalue or between the first and 
second eigenvalues. 

Eigen Solution Results Corresponding to Axial Excitation. 

Table 4 contains the first-natural-frequency eigenvector cor- 
responding to axial excitation for C r = 0.36, 0.72mm. As with 
the earlier results of Table 2, the damping factor and undamped 
natural frequency both drop as a w (0) increases. Also, the re- 
sults are relatively insensitive to changes in seal clearance. The 
undamped natural frequencies are higher for this type of mode 
shape than for the corresponding lateral-excitation mode 
shapes. 

Figure 9 illustrates the mode shapes corresponding to axial- 
excitation. Comparing Figs. 7 and 9 shows a pronounced dif- 
ference for the M,jts) eigenvectors. For the axial case, u sl (5) is 
real and, in fact, coincides with the u^s) solution. This result 
is predicted by equation (28) since u sl is uncoupled from u ei 
and p x and satisfies the zeroth-order continuity equation. The 
u S \(s) eigenvector magnitudes of Figs. 1(b) and 9(b) are similar; 
however, the magnitudes are much larger for the axial-exci- 
tation modes. The amplitude curves for /^(s) are similar in 
Figs. 7(c) and 9(c); however, the axial excitation cases are much 
larger and are relatively insensitive to changes in w w (0). 


Discussion of Approach and Results 

The eigenvalues presented in Tables 2 through 4 seem to be 
entirely consistent with the forced-response curves of Figs. 4 
and 5; specifically, the natural frequencies generally lie where 
they are expected, and they and their damping factors vary 
with Ueo( 0) as expected. The eigenvectors of Figs. 7, 8, and 9 
are difficult to comment on, given that (to the writer's knowl- 
edge) no one has tried to either calculate or measure this type 
of eigenvector in the past. The mode shapes clearly satisfy the 
boundary conditions and are consistent with the damped ei- 
genvalues to which they correspond. 

The homogeneous versions of equations (14) do not (to the 
writer) represent an obvious eigenvalue problem, and attempts 
to formulate a recognizable eigenvalue problem by conven- 
tional means were not successful. Specifically, because the 
governing equations (20) are linear, one can differentiate and 
substitute to obtain a single third-order equation in one of the 
variables. Unfortunately, this approach makes the boundary- 
condition implementation very difficult. Moreover, the final 
governing equation is itself not amenable to classical eigen- 
approaches such as central-difference finite differences; such 
approaches have proven to be notably unsuccessful for fluid 
mechanics problems. 

The approach used here to obtain eigenvalues converged 
rapidly, but the convergence characteristics suggest that the 
convergence space is not convex. Specifically, full correction 
steps predicted by the secant algorithm could not be taken. 
Generally speaking, correction steps would be reduced by a 
factor of ten at the outset. The residual error in I Del, the 
magnitude of D et would then be rapidly reduced until an ov- 
ershoot was observed. The problem would then be restarted 
at the <x t corresponding to I Del /nm, the secant correction 
factor reduced by an additional factor of ten, and a further 
rapid incremental reduction in (Del would again be realized 






until an overshoot was again experienced. By following this 
"restart” approach, the eigenvalues could be calculated to any 
accuracy desired. 

The developed algorithm would appear to be applicable to 
a range of fluid perturbation problems where the fluid me- 
chanics might have an appreciable influence on "acoustics”; 
specifically, in circumferences where convective acceleration, 
Coriolis acceleration, and/or centrifugal acceleration terms are 
not negligible. It would be adaptable to more strictly numerical 
CFD approaches such as the control-volume-based methods 
of Patankar (1980). 

The eigensolutions which were calculated in this work were, 
as specified, stimulated by forced frequency-response solu- 
tions. Obviously, more general solutions could be calculated; 
e.g., solutions of the form 

u s\ = u slc cosn6 + UjjjSin/70 
u 9 , = u eic co$n6 + u e ^in/t^ 

Pi = P\c£Osnd + p ls s\nn$ 

with n an arbitrary integer, could be assumed for eigenvalue 
solutions. The solutions developed and presented here only 
correspond to n = 0,1. Also, while the presented solutions cor- 
respond to clearance-change excitation due to impeller motion, 
these modes could be excited by strictly fluid oscillation. 

For axial excitation, Childs* (1990) work indicates that the 
impeller-shroud forces have a negligible influence on pump 
vibrations. However, Williams and Childs (1989) demonstrate 
that "fluid resonance" phenomena can have an appreciable 
influence on lateral rotordynamics of pumps. 
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Fig. 9(6) «/*,(•) flrst-natural-frequency eigenvector corresponding to ax* 
lal excitation; C,m 0.36mm 




Flfl* Pi<») flrst-natural-frequency eigenvector corresponding to ax- 
ial excitation; C r « 0.36mm 
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Influence of Impeller Shroud 
Forces on Turbopump Rotor 
Dynamics 

t-t-' 1 

\ The shrouded-impeller leakage path forces calculated by Childs ( 1987) have beer 
J analyzed to answer two questions. First , because of certain characteristics of tht 
results of Childs , the forces could not be modeled with traditional approaches . 
’ Therefore , an approach has been devised to include the forces in convent iona . 

rotordynamic analyses. The forces were approximated by traditional stiffness ; damp- 
[ ing and inertia coefficients with the addition of whirl-frequency-dependent direct 
: and cross-coupled stiffness terms . The forces were found to be well-modeled with 
. this approach. Finally , the effect these forces had on a simple rotor-bearing system 
was analyzed , and, therefore , they , in addition to seal forces , we/r applied to a 
Jeffcott rotor. The traditional methods of dynamic system analysis were modifiea 
to incorporate the impeller forces and yielded results for the eigenproblem, frequency 
response , critical speed, transient response, and an iterative technique for finding 
the frequency of free vibration as well as system stability. All results lead to the 
conclusion that the forces have little influence on natural frequency but can have 
appreciable effects on system stability. Specifically, at higher values of fuid swiri 
at the leakage path entrance, relative stability is reduced . The only unexpected 
response characteristics that occurred are attributed to the nonlinearity of the model. 


Introduction 

The vibration of centrifugal pumps has received increasing 
attention recently because of the inability of current analysis 
techniques to adequately predict the dynamic characteristics 
of pump designs. Failure to accurately predict vibrations has 
resulted in the loss of considerable amounts of money in down 
time from severe vibration problems. Massey (1985), for ex- 
ample, described an eleven-stage pump that became unstable 
when its running speed exceeded its critical speed by 25 percent. 
In other words, it whirled at 80 percent of running speed. 
Another example occurred in the High Pressure Oxygen Tur- 
bopump (HPOTP) of the Space Shuttle Main Engine (SSME) 
which also whirled at 80 percent of running speed (Childs and 
Moyer, 1985). 

The hydrodynamic forces generated in many of the fluid- 
filled gaps within the pump are well-established contributors 
to the problems cited above. This article is concerned specif- 
ically with the forces developed along the leakage path between 
the impeller and the shroud, as shown in Fig. 1. As the fluid 
is discharged from the impeller, some will return to the lower- 
pressure, suction side by way of this leakage path. Leakage is 
minimized typically with wear ring seals, as shown. Hydro- 
dynamic forces are developed along the entire leakage path, 
i.e., the shroud section as well as the seal. Note that these are 
only parts of the total force on the impeller and that the 
impeller/ volute region, balance drums, inducers, shaft, etc. 
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also contribute to the dynamics of pumps, although they are 
not considered here. 

Unfortunately, measurements of the total force on actual 
impellers have typically been obtained using pumps that have 
been modified to minimize leakage path forces. However, some 
results have been reported that are of importance here. For 
example, research at Cal Tech (Adkins, 1976) led to the con- 
clusion that the pressures in the shroud annulus contributed 
from about 50 percent to 75 percent of the total stiffness acting 
on the impeller. The pump used included separation rings and 
an enlarged shroud clearance space to minimize leakage flow 
forces. 

Bolleter et al. (1985) used a rocking-arm mechanism to ver- 
tically translate the spinning impeller. Their pump had norma! 
clearances in the leakage path; and, therefore, the forces meas- 


Wear ring seals 



Fig. 1 Typical saal configurations for a multistage pump 
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urcd were more realistic. The results demonstrated that the 
nondimensionalized cross-coupled stiffness value was about 
twice that measured on the pump at Cal Tech (Jery et al., 
1984) This suggests that the leakage path region in pumps can 
also reduce the stability of impeller motion. Bolleter et al. 
(1989) have recently presented extensive force-coefficient re- 

suits at off-design conditions. . 

Almost all of the analytical attempts at predicting leakage 
path forces have been concerned with seals. However Childs 
has extended his previous techniques in seal analysis (1980, 81 , 
82a 82b 82c) to apply them to the flow within the clearance 
Se smrounding^ shroud (1987). His work will be used 
to model the seal and shroud forces in this article and will now 

h'cSJSu^a bulk-flow approach to obtain the governing 
equations of the flow of a differential element of fluid. ln the 
seal problem, three equations were required: axial- and cir- 
cumferential-momentum equations as well as the continuity 
equation In the shroud problem, the axial-momentum equa- 
tion was replaced by a path-momentum equation, in y ° d “ c >"8 
additional terms which described the centrifugal and Coriolis 
accelerations of the fluid element. After a perturbation ex- 
pansion of the equations in the eccentricity ratio S J 
relationships were solved at various values of the wliirl fre- 
quency, fl. Integration of the resulting pressure distribution 
acting on the shaft or shroud yielded the radial and tangential 
force as functions of whirl frequency. In addition solution 
was carried out at various fluid circumferential velocities at 
the shroud entrance. This variable is of primary importance 
in system stability (cross-coupled stiffness). Figure 2 shows the 
results of the shroud problem for three different : 
inlet swirl velocity. The leakage path analyzed corresponded 
to that used by Bolleter et al. (1985) which ran at 2000 tpm. 
Figure 3 provides the dimensions of the pump and other data 
of importance. The curves in Fig. 2 define force coefficients 
corresponding to a circular orbit of amplitude and whirl fre- 
quency ratio /. The impeller force components are linearly 
proportional to the orbit amplitude A. The sharp deviations 



Fig. 2(4 Radial forca coefficient curves (Childs, 1 987) 



the radial and tangential force coefficients of Fig. 2 
used by excitation of a "centrifugal-acceleration mode nat- 
al frequency. Childs (1990) has presented results for damped 
miral frequencies and eigenvalues for these modes. 


Nomenclature 
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. kQo 
^0 


nondimensionai di- 
rect damping coeffi- 
cient 

nondimensionalized 
cross-coupled damp- 
ing coefficient 

nondimensionalized 
radial and tangential 
forces on the rotor 
(Childs, 1987)(F) 

whirl frequency ratio 

nondimensionalized 
direct stiffness coef- 
ficient 

direct stiffness coef- 
ficient which is a 
function of the whirl 
frequency ratio (F/ 
L) 

nondimensionalized 
cross-coupled stiff- 
ness coefficient 
(Childs, 1987) 
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cross-coupled stiff- 
ness coefficient 
which is a function 
of the whirl fre- 
quency ratio (F/L) 
length of seal or im- 
peller-shroud path 

(L) 


nondimensionalized 

direct mass coeffi- 
cient (Childs, 1987) 
rotor whirl radius 
(L) 

clearance at leakage 
path entrance (L) 


clearance ratio 


rotor whirl ampli- 
tude (L) 
inlet radius of 
shroud leakage path 
or seal (L) 

nondimensionalized 

tangential fluid ve- 
locity at leakage 
path or seal entrance 
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(inlet swirl ratio) 
(Childs, 1987) 
x, y = rectangular coordi- 
nates of rotor center 

(L) 

Z = complex rotor whirl 
amplitude (L) 

z=x +jy = complex rotor coor- 
dinate (L) 

0 = phase angle of whirl 
from harmonic exci- 
tation 

Q = rotor whirl fre- 
quency (1 /T) 

W|f = natural frequency of 
free vibration (1/7') 

u) s = shaft speed (1 /T) 

L = shroud leakage path 
coefficient subscript 
r = rotor coefficient 
subscript 

j - seal coefficient sub- 
script 

sh = shaft coefficient sub- 
script 

none = overall coefficient 
subscript 
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Pump data 

shaft speed = co, = 2000 rpm 
inlet clearance = Q 0 = 3.55 mm 
impeller exit angle = 22.5 0 
flowrate = 130 1/sec 
head developed = 68 m 
seven-bladed impeller 
wotting fluid is water at 80 °F 

Fig. 3 Dimension! and data for tha pump tested by Bolleter et at (Bol- 
later et al, 1985) 


In the articles cited above, Childs initially suggested the 
following conventional, linear relationship to model the forces. 
It applies for small motion of the rotor about the centered 
position. 



C c 
-c C 


(;)*[':](:) 


( 1 ) 

If circular whirl at frequency, 17, and radius, Q t is assumed 
such that 


x=Qcosfl/, (2) 

y-Q sin fl/, 

the radial and tangential force components can be stated 

After nondimensionalizing the terms as defined in Childs 
(1987), Eqs. (3) become 


j=Mf 2 -cf-K 


(4) 


where the definitions of the coefficients are included in the 
Nomenclature and the other variables are defined as follows. 

2r lap’ non< ^ mcns * on ^ zc ^ force 

_ Q 

Q * nondimensionalized clearance 

Wo 


/=-, whirl frequency ratio 

u 

F 

— , force coefficient. 

<7 

The quadratic relationships of equations (4) are to be used to 
obtain the dynamic coefficients of equations (1). Specifically, 
a least-squares curve fit of the force curves will yield the coef- 
ficients. This is an adequate approach for seal forces because 
they, in general, do follow a quadratic trend. However, the 
curves of Fig. 2 illustrate that shroud forces cannot be ap- 
proximated by equations (4); and, therefore, a traditional lin- 
ear model as expressed in equations (1) will not adequately 
represent the forces. Although these results are unexpected, 
experimental measurements made by Franz and Arndt (1986) 
of impellers with inducers yielded very similar results. 


Problem Definition 

At this point, it is not known, assuming the force curves of 
Fig. 2 are valid, how they can best be included in traditional 
rotordynamic calculations of stability, response, critical speeds, 
etc. To solve this problem, the present research has been con- 
ducted to answer the following questions. 

(/) How can the force curves found by Childs (1987) be 
incorporated in the traditional rotordynamic analyses 
of system response, stability, critical speeds, etc.? 

(2) What effects do the forces have on a rotor-bearing 
system with respect to stability, natural frequencies, 
imbalance response, etc.? 

Before answering these questions, observe that the forces pre- 
sent no new problems in a transient analysis because they could 
be applied as whirl-frequency-dependent forcing functions to 
a simple rotor model such as that developed by Jeffcott (1919). 
Although the force components are linearly proportional to 
the orbit amplitude, the equations of motion would be non- 
linear because the forces depend on the whirl frequency, de- 
fined kinematically as 

0 xy-yx 

w,(x 2 +j?y ^ 

but this could be handled by integrating the equations of mo- 
tion numerically with traditional methods. However, a tran- 
sient analysis alone is not sufficient to obtain an understanding 
of the forces, and this is why answers to the questions posed 
above are necessary. The following section describes the model 
used to analyze the shroud forces, and results are then given 
for the eigenvalue problem and frequency response, respec- 
tively. 


Rotordynamic Model 

As shown in Fig. 4, the model is a Jeffcott rotor (Jeffcott, 
1919) under the application of seal and shroud forces. The 
rotor itself represents a double-suction impeller, yielding two 
symmetrical leakage paths. In summary, the forces on the rotor 
are the shaft stiffness, two identical seal forces, and two iden- 
tical shroud forces. The mathematical representation of each 
force will now be discussed. 

A 90.7 kg (2001b) rotor is assumed. In addition, the shaft 
stiffness used (K sh * .2987 MN/m) makes the natural frequency 
of forward whirl (to be defined later) equal to 80 percent of 
shaft speed. This choice is not completely arbitrary since this 
frequency corresponds to that observed in the examples cited 
in the Introduction (Massey, 1985; Childs and Moyer, 1985). 

The seals modeled are smooth, wear ring seals with a clear- 
ance and diameter of .36 mm and 236 mm, respectively. The 
least-squares approach of Childs described earlier has been 
earned out, and the resulting dynamic coefficients are given 
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seal forces 



shroud forces 

Fig. 4 Jeffcott based, double-suction Impeller leakage path model 


Table 1 Coefficients for one wear ring seal 


VtL 

M, 

kg 

c, 

Ns/m 

c* 

Ns/m 

K a 

MN/m 

k, 

MN/m 

.5 

.3111 

3356.0 

78.6 

.6121 

.4628 

.6 

.3398 

2167.0 

89.6 

.5443 

.4611 

.7 

.3714 

970.9 

101.3 

.4755 

.4455 


in Table 1. The analysis used has at times under-predicted the 
true value of the direct stiffness. Therefore, K s is doubled when 
the coefficients of the entire model are assembled. This mod- 
ification of K s has no influence on the results, since the shaft 
stiffness is selected to yield a desired system natural frequency. 

As described in the introduction, conventional stiffness, 
damping and mass coefficients cannot completely describe the 
impeller-shroud forces calculated by Childs and shown in Fig. 
2. Therefore, the following method has been devised. First, a 
curve fit of the force curves is carried out, yielding dynamic 
coefficients which describe the qu adratic variation of the forces 
with respect to the whirl frequency ratio as described above. 
The difference between these curve fits and the actual force 
curves is modeled by direct and cross-coupled stiffness coef- 
ficients which are functions of whirl frequency ratio. The math- 
ematical equivalent of this approach will now be discussed. 

The following modified forms of equations (4) express the 
method with which the impeller-shroud forces are modeled. 

§ = Ml/ 2 ~ ?if- l* t + * * t/ )1 
Q 


7 = [£.+£•(/)], 
Q 


( 6 ) 


where R* and fc* represent the nondimensionalized deviations 
between the force curves and the approximating quadratic 
expressions. Figure 5 illustrates R * and fc * , respectively. Note 
that, although they represent forces proportional to displace- 
ment, they are not traditional stiffness coefficients since they 
are functions of whirl frequency ratio. Separating the force 
coefficients into whirl-frequency-independent and whirl-fre- 
quency-dependent terms does not introduce any approximation 
to the force definitions. Equations (1) are now rewritten as 


A K l +K'( 
[- k L -k • 


</) k L + k*if) 
( f ) Kl+K'V) 



In a mathematical sense, the tangential force accounted for 
by k* could be modeled as a damping coefficient; however, 
the result of Fig. 5 (b) confirms the wisdom of modeling the 
deviation of the tangential-force coefficient by a whirl-fre- 
quency-dependent cross-coupled-stiffness-coefficient. Specif- 
ically, Jt* oscillates about zero and approaches zero as the inlet 
swirl ratio is reduced. Finally, observe that equations (7) are, 




Table 2 Approximating coefficients for one leakage path 



C L 

CL 

Kl 

ki 

kg 

Ns/m 

Ns/m 

MN/m 

MN/m 

3.469 

2037 

969.9 

.06781 

.3137 


in general, nonlinear differential equations. To understand this 
statement, recall that the general definition of the whirl fre- 
quency ratio,/, is given in equation (5). If the rotor precession 
has a constant radius and rate, / is a constant and the motion 
is linear. However, for general motion, / is variable and a 
function of the motion. 

After performing the asymptotic least-squares curve fit to 
Fig. 2, the resulting coefficients of equation (6) are given in 
Table 2. 

In summary, the following equations represent how the over- 
all coefficients for the model are defined. 

A/=A/, + 2(A/, + A/ l ) 

C=2(C, + C £ ) + C 

c=2(c, + c £ ) (8) 

K=K sh + 2i2K s + KiJ + 2K'=K + 2K' 
k = 2(k s + k L ) +2A:* =k + 2k*. 

Particular coefficients are doubled because there are two 
seals and mo leakage paths in a double suction pump. The 
quantity, C, must be added to the model to yield reasonable 
stability. If the eigenvalues of the overall system are obtained 
without adding C, the system is found to be unstable. There- 
fore, by adding another damping term (C- 10,422 Ns/m), 
reasonable stability, which is consistent with operating pump 
experience, is ensured. The following equation represents the 
complete model used to model the forces on the rotor. 
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Table 3 Coefficients for overall model 


Table 4 Natural frequencies of free vibration for the pump model with 
and without K * and k m 


UtL 

M 

C 

c 


mum 


kg 

Ns/m 

Ns/m 



.5 

98.240 

21208 

2097.1 

2.8827 

1.5549 

.6 

98.298 

18830 

2119.1 

2.6115 

1.5494 

.7 



2142.4 

2.3363 

1.5184 



Fig. 6 Real part of the eigenvalues for the forward-whirl mode 



Whirl Frequency Ratio (ND) 

Fig. 7 Imaginary part of the eigenvalues for the forward-whirl mode 


M 0 
0 M 




K{f) k(f) 
k(f) K(f ) 



Table 3 contains the numerical values for this analysis. These 
coefficients combine to make the damping for the uncoupled 
system, defined by C/2\jKM t equal to 63 percent, 59 percent, 
and 54 percent of critical damping for the three swirl ratios 
shown, respectively. Again, R* and k* are illustrated in Fig. 
5 and are defined analytically with cubic-spline curve fits. 


Natural Frequency Determination 

The characteristic equation for Eq. (9) is 

[Ms 1 + Cs + K(f)l?+[cs+k(f)\ 1 = 0, (10a) 

where 5 = $/oj,. For each whirl-frequency ratio/, two pairs of 
complex conjugate roots of the form 

s = Of i— 1,2 (1 Ob) 

can be calculated corresponding to forward and backward 
whirling modes. The concern here is finding the natural fre- 


UtL 

— with 
K* and Jfc* 

— without 

w, 

K * and k * 

Percent 

Decrease 

.5 

.747 

.727 

2.68% 

.6 

.777 

.724 

6.82% 

.7 

.800 

.715 

10.63% 


Tabie 5 Real parts of the eigenvalues for the model with and without 
K m and k * 


UtL 

Real part with 
K * and k* 
1/sec 

Real part without 
K * and k * 
1/sec 

Percent 

Decrease 

.5 

-55.9 

-60.3 

7.87% 

.6 

-31.1 

-47.2 

51.80% 

.7 

-13.8 

-34.5 

150.00% 


quency for forward whirling motion. Figure 7 illustrates the 
eigenvalue solutions to Eq. (10a) versus/. 

Free whirl of the rotor occurs only when the imaginary part 
of an eigenvalue equals the whirl frequency ratio /in Fig. 7; 
i.e., X 2 =/. The three points at which the X 2 =/Iine intersects 
the X(/) curves defines the forward-whirl natural frequencies 
for the three inlet swirl ratios. The real part of the eigenvalue 
at this natural frequency defines the rotor’s damping. 

To find the “natural frequency” and damping of free whirl, 
the following iterative technique has been used. First, a whirl 
frequency ratio is assumed. From this value, the corresponding 
values of K and k are obtained from the data shown in Fig. 
5 from which the eigenvalues are found using Eq. (10). These 
first two steps are equivalent to finding a point on the curves 
of Figs. 6 and 7. The imaginary part of the forward-whirl 
eigenvalue is compared to the whirl frequency used. If they 
are the same, the natural frequency at which free vibration 
takes place is defined. If they are different, the imaginary part 
of the forward-whirl root becomes the assumed whirl frequency 
ratio, and the same steps are carried out until convergence 
occurs. Note that this calculation procedure yields the damped 
natural frequency for the system. 

Tables 4 and 5 show results of the above algorithm. Also 
included are the eigenvalues of the same model neglecting K * 
and k •, which correspond to a pure quadratic approximation 
of the leakage path forces as in equations (4). The eigenvalues 
illustrate at least two important conclusions about the effects 
of the values of K* and k*. First, the natural frequency in- 
creases when K* and k • are included as well as when the swirl 
ratio is increased. The variation of K * in Fig. 5 explains these 
observations. More importantly, the percent differences be- 
tween the two models show that the values of K* and k * have 
very little effect on natural frequency. Second, stability is re- 
duced when K* and k • are included and when the swirl ratio 
is increased. The variation of k • in Fig. 5 is the cause of these 
observations. Finally, the percent differences between the two 
models shows that K* and k * can cause appreciable reductions 
in relative stability. 

To show that the stability and frequency of free vibration 
are determined only by the root obtained from the iterative 
approach described earlier, the equations of motion (9) were 
integrated in a series of transient, free-vibration simulations. 
Initial conditions were an initial displacement of .127 mm (.005 
in) and velocity equal to Riw s (synchronous whirl). Figure 8 
includes the resulting orbit and a time history of the whirl 
frequency ratio for free vibration and an inlet swirl ratio of 
0.7. As shown, the rotor executed a well-damped spiral orbit. 
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F)g. 8(«) Rotor orbit during free vibration 



Fig. 8(b) Time history of rotor whirl frequency ratio during free vibration 



WTuH Frequency Rino (ND) 

Fig. 9 Real part of the eigenvalues for the forwsrd'Whirt mode with a 
12.5 percent reduction in direct damping 


The whir! frequency ratio of the orbit was exactly 0.8, which 
corresponds to the natural frequency found earlier. The other 
two swirl ratios yield similar results. 

To verify that the stability of the mode! is described only 
by the real part of the eigenvalues obtained with the iterative 
approach described above, the direct damping coefficient was 
reduced by 12.5 percent producing a range of unstable eigen- 
values as shown in Fig. 9 where the real part of the forward- 
whirl root is positive between whirl frequency ratios of 0.41 
and 0.57. The natu ral freque ncy of this altered model c orre- 
sponds to a stable eigenvalue and is .825 w,. Integration of the 
eq uations o f motion wer e altered to include a harmonic forcing 
function with an amplitude of 5N which forced the rotor in a 
circular direction at a frequency of .5w,. Figure 10 includes 
the resulting orbit, whirl frequency ratio time history, and a 
Fast Fourier Transform (FFT) of the ^-coordinate, respec- 
tively. Even though the forcing frequency was within the un- 
stable range of frequencies, the figures show that the system 
was stable in that its orbit did not grow without bound. The 
loops and limit cycle behavior are the results of the nonlin- 
earities in the model. The FFT shows that the steady state orbit 
consisted of response at both the forcing frequency, ,5w Si and 
the natural frequency, .825 w s . Other forcing frequencies were 
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Fig. 10(a) Rotor orbit during harmonic excitation at 50 percent of shaft 
speed and 12.5 percent reduction in direct damping 



Fig. 10(b) Time history of rotor whirl frequency ratio during harmonic 
excitation at 50 percent of shaft speed and a 12.5 percent reduction in 
direct damping 



Fig. 10(c) FFT of the response of the rotor during harmonic excitation 
at 50 percent of shaft speed and a 12.5 percent reduction in direct 
damping A 


also tried. Using a stable value of 0.7 w st for example, the 
response did not behave in the same manner; instead, the orbit 
reached a steady state circular path at the excuatTon frequency 
as would be expected from a linear system. Howe ver, forcing 
at .315w st the characteristics of the results from forcing at ,5w s 
occurred again, as can be seen in Fig. 11, which ag ain contain 
the orbit, tune history of the whirl frequency ratio and an FFT 
of the response. Only forcing frequencies below and within 
the unstable zone resulted in this behavior. Only when the 
direct damping coefficient was reduced until the natural fre- 
quency corresponded to an unstable eigenvalue did the re- 
sponse become unstable in a linear sense, where the orbit grew 
without bound at the natural frequency. The fundamental 
result from this analysis is a verification that a range of whirl 
frequencies yielding eigenvalues with positive real parts does 
not cause the system response to grow without bound unless 
it includes the system natural frequency. In addition, the non- 
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Fig. 1 1<j) Rotor orbit during harmonic excitation at 37.5 percent of shaft 
speed and a 12.5 percent reduction In direct damping 
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Fig. 12 Frequency response: whirl amplitude 



Fig. 11(b) Time history of rotor whlri frequency ratio during harmonic 
excltallon at 37.5 percent of shaft speed and a 12.5 percent reduction 
In direct damping 



Whirl Frequency Rauo (ND) 



Fig. 11(c) FFT of the response of the rotor during harmonic excitation 
at 37.5 percent of shaft speed and a 12.5 percent reduction in direct 
damping 


linearities of the system can cause a limit cycle to occur under 
certain conditions. 

Finally, the critical speed of the rotor cannot be calculated 
unless force curves are available at other shaft speeds. In this 
case, the natural frequency is calculated with the iterative 
method just described, but at a number of shaft speeds. A 
plot of natural frequency as a function of shaft speed is then 
drawn, analogous to “critical speed maps’* in fluid-film bear- 
ing analysis. The critical speed is defined as the speed where 
the natural frequency is equal to the shaft speed. 


Frequency Response 


As discussed earlier, the following are the equations of mo- 
tion of the model 




Fig. 13 Frequency response: phase angle 


where the right-hand side represents a harmonic forcing func- 
tion appropriate to frequency response analyses. As in tradi- 
tional methods, the transient solution is assumed to decay with 
time, and only the steady-state response of the rotor is con- 
sidered. Obviously, frequency response solutions are only valid 
when the system is stable, precluding the type of response 
illustrated in Figs. 10 and 11. The forcing frequency is defined 
by Q which has been used to represent whirl speed. Since the 
steady-state motion is assumed to occur at the same frequency 
as the excitation, this should not cause confusion. Note that 
the model is now linear because the whirl frequency ratio and, 
hence, the stiffnesses are constant. Multiplying the second 
equation by the complex number /, the equations can be com- 
bined into one in terms of the complex variable z. 

Mz+(C-ic)z + (K-ik)z = Fe /ai . (12) 

Assuming steady-state motion of the following form 

z = Re?+, (13) 

where R denotes the whirl amplitude, and <t> is the phase angle 
yields 


R K 

F/K Vl*- + cQ ] 2 + [cn - Ar] J (14) 

, . _.f *-CP 1 
* 811 [at-a«i 2 +cqJ' 

Equations (14) were solved using the present model data, 
and the results for a swirl ratio of 0.7 are shown in Figs. 12 
and 13. Also included are the same results for the model with 
K* and k • neglected. Again, in obtaining frequency curves as 
shown in Figs. 12 and 13, the stiffness terms must be calculated 
at the corresponding whirl frequency ratio. 

The swirl ratio of 0.7, used in Figs. 12 and 13, exhibited the 
largest variation in amplitude and phase at the whirl frequency 
ratio of .80, the system natural frequency. An amplitude dif- 






fcrcncc of about 60 percent is shown between the curves, sug- 
gesting that K* and k * have significant influence on the system 
primarily at this frequency. The swirl ratios of 0.6 and 0.5 
exhibited similar but decreasing results in both the amplitude 
and phase curves, having differences in amplitude at the peak 
value of 35 percent and 6 percent, respectively. In fact, the 
curves were almost identical in the case corresponding to a 
swirl ratio of 0.5. Note that these resonance peaks occurred 
in all cases at the corresponding natural frequencies found in 
the previous section, further verifying the preceding results. 
The difference in the response to synchronous imbalance ex- 
citation (/= 1) between the two models is small, and this sug- 
gests that the leakage path forces have negligible effects on 
imbalance response. However, the presence of a subsynchron- 
ous excitation near the natural frequency can result in large 
increases in response amplitude when k* is included. 

Summary and Conclusions 

The impeller-shroud forces have been separated into (a) con- 
stant (whirl frequency-independent) stiffness, damping and 
mass coefficients and ( b ) direct and cross-coupled added stiff- 
ness coefficients which are functions of whirl frequency. With 
this approach, the model can be analyzed using traditional 
techniques with a few modifications. A new iterative technique 
was used in the determination of the rotor’s natural frequency 
and damping. The whirl-frequency-dependent stiffness terms 
were seen to be convenient and effective models for the im- 
peller-shroud force nonlinearities. 

When the values of the added stiffness coefficients, a and 
£•, were included in addition to the frequency-independent 
coefficients, a number of effects resulted. First, the natural 
frequency was seen to increase by a small amount, suggesting 
that K * could be neglected when only the natural frequency 
is important. In addition, the stability is seen to decrease when 
k* is included. This effect is appreciable, so stability calcu- 
lations should include k *. The frequency response analysis 
showed synchronous response to imbalance is increased by the 
stiffnesses although by only a small amount. Appreciable dif- 
ferences in response due to subsynchronous excitation occur 
only at or near the natural frequency and become negligible 
at low values of the inlet swirl ratio. 

The nonlinearity introduced by the whirl-frequency-de- 


pendency of the force coefficients resulted in limit cycles and 
nonsynchronous response in some cases. 


References 

Adtkins, D. ? 1976, “Analyses of Hydrodynamic Forces on Centrifugal Pump 
Impellers ” M.S. Thesis. California Institute of Technology, Pasadena. 

Bolleter, U., Wyss, A., Write, I., and Sturchler, R., 1985, “Measurements 
of Hydrodynamic Interaction Matrices of Boiler Feed Pump Impellers,” ASME, 
Paper No. 85-DET-I48. 

Bolleter, U., Leibundgut, E., Sturchler, R., and McGoskey. T., 1989, “Hy- 
draulic Interaction and Excitation Forces of High Head Pump Impellers," pre- 
sented at the ASCE/ASME Pumping Machinery Symposium, July 9-12, 1989. 

Childs, D., Dressman, J., and Childs, S., 1980, “Testing of Turbulent Seals 
for Rotordynamic Coefficients,” Proc . of Workshop on Rotordynamic Insta- 
bility Problems in High Performance Turbomachinery, Texas A&M University, 

NASA CP-2133, pp. 121-38. , , » 

Childs D 1981, “Convcrgcnt-Tapered Annular Seals: Analysis for Rotor- 
dynamic Coefficients,” Symposium Vol., Fluid /Structure Interactions Turbo- 
mach., ASMW Winter Annual Meeting, pp. 35-44. 

Childs, D., 1982a, “Dynamic Analysis of Turbulent Annular Seals Based on 
Hirs’ Lubrication Equations,” ASME Paper No. 82-Lub-4l. 

Childs, D., 1982b, “Finite-Length Solutions for Rotordynamic Coefficients 
of Turbulent Annular Seals,” ASME Paper No. 82-Lub-42. 

Quids, D., and Dressman, J., 1982c, “Testing of Turbulent Seals for Ro- 
tordynamic Coefficients,” Proceedings of Workshop on Rotordynamic Insta- 
bility Problems in High Performance Turbomachinery , Texas A&M University, 
NASA CP-2250, pp. 157-71. 

Childs, D., and Moyer, D., 1985, “Vibration Characteristics of the HPOTP 
(High Pressure Oxygen Turbopump) of the SSME (Space Shuttle Main Engine),” 
ASME Journal of Engineering for Gas Turbines and Power, Vol. 107, pp. 152- 


159. 

Childs, D., 1987, “Fluid-Structure Interaction Forces At Pump-Impeller- 
Shroud Surfaces for Rotordynamic Calculations,” Texas AAM University, Me- 
chanical Engineering Dept., ASME Vibration Conference, Boston, Mass. 

Childs, D., 1990, “Centrifugal-Acceleration Modes for Incompressible Fluid 
in the Leakage Analysis Between a Shrouded Pump Impeller and Its Housing," 
ASME Journal op Vibration and acoustics, May 1990. 

Franz R and Arndt, N., 1986, “Measurement of Hydrodynamic Forces on 
the Impeller* of the SSME,” Report No. E249.2, California Institute of Tech- 
nology, Pasadena. 

Jeffcott H 1919, “The Lateral Vibration of Loaded Shafts in the Neigh- 
borhood of a Whirling Speed-The Effect of Want of Balance,” The Philo- 
sophical Magazine, Vol. 6, No. 37, pp. 304-314. 

Jery, B., Acosta. J.. Brennen, C., and Caughey. T., 1984, “Hydrodynamic 
Impeller Stiffness, Damping and Inertia In the Rotordynamics of Centrifugal 
Flow Pumps,” Proceedings of Workshop on Rotordynamic Instability Problems 

rr> l n T#vai AAM UnivCfSitV. NASA CP-2338, 


Massey. I. C.. 1985, “Subsynchronous Vibration Problems in High Speed 
Multistage Centrifugal Pumps,” Proceedings of the 14th Turbomachinery Sym- 
posium, Turbomachinery Laboratories, Texas AAM University, pp. 11-16. 


Journal of Vibration and Acoustics 


OCTOBER 1991, Vol. 113/515 


^ n 

»— » :r 
vO H- 
vO H*‘ 

fO Cl. 

cr CO 

<• 




i 

i 


feifi 




Reprinted from March 1992, Vol. 114, Journal of Fluids Engineering 




D. W. Childs 

Turbomachinery Laboratories, 
Mechanical Engineering Department, 
Texas A&M University, 
College Station, Texas 77843 


Pressure Oscillation in the 
Leakage Annulus Between a 
Shrouded Impeller and Its 
Housing Due to Impeller- 
Discharge-Pressure Disturbances 

An analysis is presented for the perturbed flow in the leakage path between a 
shrouded-pump impeller and its housing caused by oscillations in the impeller- 
discharge pressure. A bulk-flow model is used for the analysis consisting of the path- 
momentum, circumferential-momentum, and continuing equations. Shear stress at 
the impeller and housing surfaces are modeled according to Hirs’ turbulent lubri- 
cation model. In the present analysis, perturbations of the impeller discharge pressure 
are used to excite the fluid annulus. The circumferential variation of the discharge 
pressure is expanded in a Fourier series up to order n where n, is the number of 
impeller blades. A precession of the impeller wave pattern in the same direction or 
opposite to pump rotation is then assumed to completely define the disturbance 
excitation. Predictions show that the first ( lowest-frequency ) “centrifugal-accel- 
eration" mode of the fluid within the annulus has its peak pressure amplitude near 
the wearing-ring seal. Pressure oscillations from the impeller can either be attenuated 
or (sharply) magnified depending on: (a) the tangential velocity ratio of the fluid 
entering the seal, ( b ) the order of thefourier coefficient, and (c) the closeness of 
the precessional frequency of the rotating pressure field to the first natural frequency 
of the fluid annulus, and (d) the clearance in the wearing-ring seal. 


Introduction 

The present work is stimulated by experiences with the SSME 
HPFTP (Space Shuttle Main Engine, High Pressure Fuel Tur- 
Bopump) wearing-ring seaTs. A stepped, 3-cavity, tooth-on- 
rotor, labyrinth-seal design is used. The stator for the seal is 
made from KEL-F, a plastic' that is somewhat similar to nylon. 
In some cases, post-test inspection of the stator element has 
revealed that interior points in the stator material have melted 
and then resolidified, despite being in contact with liquid hy- 
drogen. One hypothesis for this exceptional outcome was that 
the material had been subjected to cyclical stresses which gen- 
erated heat due to hysteresis. Because of poor conduction 
properties of the material, the heat could not be dissipated, 
the temperature rose, and melting resulted. “What pressure 
oscillations are driving the cyclical stresses?”, is an obvious 
question in reviewing this scenario. The present analysis ex- 
amines “centrifugal-acceleration” modes, arising between the 
impeller shroud and its housing and driven by pressure oscil- 
lations from the pump, as an answer to this question. 


‘The work reported herein was supported by NASA Marshall Space Right 
Center under contract NAS 8-37821 ; contract technical monitor: James Cannon. 

Contributed by the Fluids Engineering Division for publication in the Journal 
of Fluids Engineering. Manuscript received by the Fluids Engineering Division 
August 20, 1990. 
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Figure 1 illustrates an impeller stage of a multi-stage cen- 
trifugal pump. Leakage along the front side of the impeller, 
from impeller discharge to inlet, is restricted by a wearing-ring 
seal, while leakage along the back side is restricted by either 
an interstage seal or a balance-piston discharge seal. The pres- 
ent analysis considers perturbed flow in the leakage paths be- 
tween the impeller-shroud surface and its housing. 

Prior analyses by the author of those annulli have been 
concerned with lateral (1987, 1989) and axial (1990a) reaction 
forces developed by the impeller shrouds as a consequence of 
harmonic clearance changes due to impeller motion. These 
analyses have been based on “bulk-flow” models which neglect 
the variation in the dependent variables across the fluid film. 
The model consists of the path and circumferential momentum 
equations and the continuity equations. 

The analyses cited Tiave yielded force and moment coeffi- 
cients due to impeller motion but have also predicted “reso- 
nance” phenomena, which are caused by the centrifugal- 
acceleration body forces present in the path momentum equa- 
tions. An algorithm was developed (1990b) to calculate the 
complex eigenvalues and eigenvectors associated with these 
resonances. In the present analysis, the harmonic response of 
the flow within the annulus is examined due to time and cir- 
cumferential variations in the discharge pressure of the im- 
peller. 
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element of thickness /f(S, 0, /). From Childs (1987), the con- 
tinuity equation can be stated 


dH d 


1 0 


H\ dR 


li + Ts^ H)+ R3e lu ^ ) + [j)Ts u ^ 0 (3a) 

where U 5 and U e are the path and circumferential bulk-velocity 
components, respectively. The path and circumferential mo- 
mentum equations are stated 

ap „u]dR 

- H Ts =T “ +T ”- pH TTs 


Jau, au,u, ju, ' 


m 


H dP 

R dd~ Te, + Ter 


Following Hirs’ approach (1973), the wall shear-stress defi- 
nitions in these equations can be stated 


ns 


= J pU]R?[\+(U t /U,) 2 } 2 


T;r=- pU]R7 r \ 1 + [{U,-Ru)/U S Y\ 2 


ms+ 1 


Tes = ™ pUeU^ll + (U e /U s f\ 2 


r», = - pU s (U»-Rw)RT[ 1 + l(U,-Rw)/U s ] 2 \ 


mr + 1 
2 


(4) 


where 


R s =2HU s /y 


(5) 


Nondimensionalization and Perturbation Analysis 

The governing equations define the bulk-flow velocity com- 
ponents (C/j, U 9 ) and the pressure P as a function of the co- 
ordinates (R6 y S) and time, /. They are conveniently 
nondimensionalized by introducing the following variables 

Us=Us/V h Ue^Ue/R^y p = P/f>Vj 

s = S/L Si r = R/R t y b = V { /R,w 

r = cj/, T= L s /V, (6) 

The present analysis examines the changes in (u st u 9y p) due 
to changes in the impeller’s discharge or inlet pressure. Fol- 
lowing conventional notation, pressure drops at the annulus 
inlet and exit are stated 

p s (o,t) - P(o,e,t) = P ( 1 + i)u)( 0 ,e,t )/2 

P(L S , 6, 0-P t (6j) =pC dt U 2 (L s ,d,t)/ 2 (7) 

Note specifically that the (upstream) supply and (downstream) 
exit pressure are now functions of time. Assume that the os- 
cillations consist of a small perturbation of the form 

PAW + *P S 1 ($,t)J>AW + *Pel (W 

and introduce nondimensional variables to yield the following 
zeroth 

A>(0)=Ao-(l-*)/2 

#*1)-A.+ C*i&(1)/2 (8) 

and first-order equations 

Ps\ (0,f ) ~P\ (0,0,/) = ( 1 + £)w„ (0,0,/) 

Pi ( 1 At) -Pei (0,/) = cw l)u„ ( 1 ,0,/) (9) 



The perturbed supply and exit pressures p s , (0,/), p e i(0,t) can 
now be specified functions of time and provide excitation for 
the perturbed flowfield within the annulus. 

Expansion of the dependent variables of Eq. (3) in pertur- 
bation equations yields: 


Zeroth-Order Equations 

(a) Path-Momentum Equation 


dp 0 1 

/dr\ 


(a r + a\ 

J_ dho J dr 

ds r 

[ds) 

v + [ 

\ * ) 

h 0 ds r ds 


( b ) Circumferential-Momentum Equation 


^ £ + ^ 2 ^ + I <rr +< W»o] /2 = 0 (10*) 

(c) Continuity Equation 

rh Q u so = 1 (10c) 

The quantities a s and a r are defined by 

o s = (L s /H 0 )\ S y <j r =(L s /H 0 )\ r (11) 

where \ s and \ r are the dimensionless stator and rotor friction 
factors: 


ms * 1 

K = nsR£U + (i u to /bu so ) 2 ]~ 

mr + 1 

K = nr/C 1 1 + [(u So - r) /bu so ] 2 ) ~T~ 

The continuity equation has been used to eliminate du so /ds 
from Eq. (10a). The momentum equations define the pressure 
and velocity distributions fora centered impeller position. They 
are coupled and nonlinear and must be solved iteratively. The 
initial condition for u 9o ( 0) is obtained from the exit flow con- 
dition of the impeller. Zeroth-order pressure boundary con- 
ditions are provided by Eq. (8). 

Figure 3 illustrates the pump-impeller and shroud geometry 
used by Bolleter et al. (1987) in their test program for radial 
force coefficients. Their tests were at best efficiency point 
(BEP) with the pump running at 2000 rpm, while developing 
68m of head and 1301/s of flow rate. The impeller has seven 
blades and an impeller exit angle of 22.5*. The test fluid is 
water at 26.6*C. For the present study, A P across the impeller 
is assumed to be 70 percent of the total head rise of the stage. 
Based on pitot-tube measurements, impeller-exit-tangential ve- 
locity is about 50 percent of the impeller discharge surface 
velocity; hence, h 0.5. 
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Table 1 Zeroth-order-solution results; C { 


0.36mm 


u 9o (0) 

Rso 

C de 

m (kg/sec) 


0.5 0.6 

9377 8907 

2.068 2.098 

4.448 4.225 


3.5mm, C r = 


0.7 

8426 

2.130 

3.997 


Table 1 provides the zeroth-order solutions for new (original- 
clearance) exit wearing-ring seals. The R& values on the order 
of 10,000 are low in comparison to the circumferential Reyn- 
olds number R m = 2 HU m /v which varies along the path, but 
is on the order of 250,000. These leakage results were obtained 
iteratively starting with guessed values for the seal inlet values 
of pressure and tangential velocities and yielding an initial 
estimate for C^. With this estimate for C*, the leakage is 
calculated through the annulus which yields new inlet condi- 
tions for the seal. Solutions are “bounced” back and forth 
between the seal and the annulus until the same leakage value 
is obtained for both flow paths (four-place accuracy). 


First-Order Equations 

(a) Path-Momentum Equation 


= 0 (12a) 


^r + u e iA2s + u 5l A) 5 

ds 

F -t u *o dUsi dus_ 1 

[dr r 36 ds 

( b ) Circumferential-Momentum Equation 

L L s 1 dp x . 

b -r rr + Ue\Aie + u s \A^ 

R { r ad 

[ -du e i 

wr -r“ + «/ — rr + u so -r 

3t r 36 os 

(c) Continuity Equation 

3u s[ oj T du 0 j (\ dr 1 dh 0 \ 

New coefficients in these equations are defined in, Childs 
(1987). 


= 0 (12 b) 


(12c) 


u s \ = u 5lc +ju slst u e] = Ue\c +ju e \s P\=Pic+jP\s U 6 ) 

reduces these real equations to three, complex, ordinary dif- 
ferential equations 


d_ 

ds 


r u > 
ue i 

A 


+ [A(nJ f s)]\u 9] =0 


(17) 


where 




B - jnmt/r 

Am/ u so (Aie+jTT)/UsQ 

l nuT 

A is - Buyo+jTT [Ais+J —— Mjo 


-J 


0 

n b l Ls 
ru s o \/W 
0 


1 rfr 1 dh 0 (fm ^ \ 

r ds h 0 ds \ r / 


(18) 

(19) 


Since there is no right-hand side to Eq. (17), the homogeneous 
solution is the complete solution and can be stated as follows 
in terms of the transition matrix and initial conditions 


( 20 ) 


The inlet initial condition u ei ( 0) is set equal to zero, and 
calculation of ^,(0) and p x ( 0) in terms of the specified bound- 
ary conditions is the immediate problem at hand. Substitution 
from Eqs. (13) and (14) into Eq. (9) yields 

n . _ n.m\ = n 4- 


r u s[ y 


u$ 1 

= [$( «,/,*)] |m„(0)| 


\ {p i(0) J 


where 

Ps\ ” Pslc +JPs U 

Pe\ =Pt\c+JPe\s (22) 

From Eq. (20) 

«,i(l) = *u(l)w,i(0) + *u(\)Pi(Q) 


Solution Procedure: First-Order Equations 

The functions p sl (6, t ), p eX {d : /) provide the boundary ex- 
citation for the first-order equations. The general form for the 
excitation takes the form 

p sl (d,t)=e jW {p s]c cos nd+p 5 i s sin nd) 

p t \ (0*/) ® e jUt (Pe\c cos n0+Peis sin nd) (13) 

where n can reasonably be expected to vary from zero (plane 
wave) upwards through multiples of the number of blades in 
the impeller. The form of Eq. (13) suggests that the 6 variation 
in boundary pressures is defined in an impeller-fixed coordinate 
system, which is precessing at the frequency Q. 

The 6 and time dependency of the dependent variables is 
eliminated by assuming the comparable, separation-of-vari- 
able, solution format 

u sX = e J/T (u s i c cos nd + w jl5 sin nd) 

u $x =e J/T (u 9 ic cos nd + u 0Xs sin nd) 

p x =e jJr (P\ c cos nd+p\ 5 sin nd) (14) 

where the coefficients are solely functions of 5, and 

/= fl/w (15) 

is the normalized precession frequency. Substituting into Eqs. 
(12) and equating like coefficients of cos nd and sin nd yields 
six first-order equations in 5. Introducing the complex variables 


/7j(l) = 4*3 j (l)Wji(O) + *330 )Pi(0) (23) 

Hence, from Eq. (21) one obtains 




Inversion of this equation yields 

Cwn(0)l _ lAl Z\2 

U>i(0)j lAl ^22. 

Jp sX Sind p eX cannot be specified independently, and a relation- 
ship between the two cannot be established without a knowl- 
edge of the fluid system beyond the current terminating orifices. 
For the purposes of this discussion, the arbitrary choice 

Psl = U Pe 1=0 

is made to examine the influence of pressure perturbations at 
the impeller exit (annulus inlet). The resulting set of initial 
conditions for Eq. (20) is then 

&£)-£]• *«-• 

The complete solution along the impeller is found by evaluating 
Eq. (20) for s€[0, 1]. 
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Numerical Results. Childs’ ( 1 990b) analysis yields complex 
eigenvalues 5 and eigenvectors u s \(s)> u$\(s) y p\(s) for the 
system modeled by Eq. (12), Figure 4 illustrates the amplitude 
and phase for the first (lowest-natural-frequency) eigenvalue 
for n = 1; « w (0) = 0.5, 0.6, and 0.7. Observe that the peak- 
pressure amplitudes lie near the wearing-ring seal, which is 
consistent with the internally-melted HPFTP wearing-ring-seal 
results cited earlier. Forced, harmonic-response solutions de- 
veloped in this study due to impeller-discharge-pressure oscil- 
lations also show the largest pressure oscillations to occur near 
the wearing-ring seal (5 = 0.95); hence, results presented here 
focus on p t (0.95). This peak-pressure-oscillation location is 
very near the exit-wearing-ring seal of Fig. 3. The first question 
to be addressed here is, “How does p\( 0.95) depend on n,/, 
and C,?” 

Figure 5 illustrates the amplitude and phase of pj(0.95) versus 
/for C r - 0.36mm, n = 5, and = 0.5. As will be explained 
later, the choice n = 5, arises because of the number of impeller 
and diffuser blades used in Bolleter et al.’s (1987) pump. The 
phase results indicate that numerous resonances exist for po- 
sitive values of /. However, only the first resonance experiences 
significant amplification. The remaining fluid modes are heav- 
ily damped. 

Figure 6 illustrates Ip, (0.95) I for C r = 0.36mm, n = 5, and 
w#o(0) = 0.5, 0.6, and 0.7. The peak-response frequency in- 
creases as Um( 0) is increased from 0,5 to 0.6, and a secondary 
peak appears around / = 2.5. Increasing w w (0) from 0.6 to 
0.7 causes an additional peak to appear. 

Figure 7 illustrates Ip, (0.95) I for C r = 0.36mm, u^O) = 
0.5, and n - 0, 1, 3, 5, and 7. The response is heavily damped 
for n = 0, rises sharply as n is increased to one, but then 
remains relatively constant as n ranges upwards over 3,5, and 
7. Figure 8 repeats the results of Fig. 7, except for worn clear- 
ances; i.e., C r = 0.72mm. Comparisons of Figs. 7 and 8 show 
that doubling the clearances reduces pressure amplification and 
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Fig. 5 (a) lft(0.95)l versus /; C r * 0.36mm, *W0) = 0.5, n = 5 
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Fig. 5(b) Phase £,(0.95) versus f ; C r * 0.36mm, u*(0) = 0.5, n m 5 
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Fig. 6 I ft (0.95) I versus f;C,m 0.36mm, n = 5, and utn(0) = 0.5, 0.6, 
and 0.7 

slightly elevates the peak-amplitude-excitation frequency/’. 
Table 2 shows f* versus n for the new(C r = 0.36mm) and 
worn seals(C r = 0.72mm). Note that/* increases more-or-less 
linearly with increasing n. 

The questions which now arise are: In a real pump, wbat 
value of n is likely to arise in impeller-pressure-discharge pat- 
terns, and what precession frequency is most likely to be present 
and dominant? Answers to these questions have been provided 
by Bolleter (1988), who presents an analysis for the pressure 
waves developed by the interaction of impeller and diffuser 
vanes or impeller vanes and volutes. For an impeller with rt } 
vanes and a diffuser with /I 2 vanes, Bolleter shows that a ro- 
tating pressure wave is developed around the impeller exit with 
n - l«j - n 2 \ diametral nodes. If #ij > n 2 the pressure wave 
rotates in the direction of the pump with the frequency n\u/ 
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Table 2 / * (peak-excitation-amplitude frequency) versus n 
for w*>(0) = 0*5 


n 

1 

2 

3 

4 

5 

6 

7 


0.3 

0.8 

1.4 

2.1 

2.9 

3.6 

4.2 

C r = 0.72mm 


0.9 

1.6 

2.4 

3.1 

3.8 

4.4 
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Fig. 7 lpi(0.95)l versus /for n * 0, 1, 3 t 5, 7; C, » 0.36mm, and u«(0) 

« 0.5 
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Fig. 8 [£,{0.95)1 versus / for n = 0, 1, 3, 5, 7; C, * 0.72mm, and u^0) 

* 0.5 

n. If n x < n 2 the precessional frequency is Note that 

n t the number of diametral nodes, cited by Bolleter is the same 
n used in Eq. (13) for the pressure excitation. Further, in terms 
of Eq. (13), ft = i/ijuj/lfl! - n 2 1 . Tyler and Sofrin (1962) 
earlier developed this same result in analyzing the noise gen- 
erated by the interaction of a rotor and stator in axial com- 
pressors of gas turbines. 

Bolleter et al.’s pump (1987) used 12 diffuser blades, and 
the impeller of Fig. 3 has 7 blades. Hence, from Bolleter (1988), 
n - 17 - 121 = 5, and ft = -7w/5 = - 1.4gj. If the impeller 
were mounted in a double volute, n = 17-21 =5, and ft 
= 7a>/5 = l,4<i>. From Fig. 5 (n = 5, u w (0) = 0.5), ampli- 
fication for / = - 1.4 and 1.4 is 0.23 and 0.5, respectively. 
Hence, pressure disturbances from the impeller would generate 
pressure oscillations about twice as large in a double volute as 
in a 12-vaned diffuser. However, in either case, because the 
predominant frequency is well removed from the peak-am- 
plitude-excitation frequency f* = 2.9, impeller pressure dis- 
turbances would actually be attenuated by the annulus. 

From BoIIeter’s equations, and the results of Figs, 5 through 
7, significant amplification of impeller-discharge-pressure var- 
iations will only arise when the number of impeller blades 
exceeds the number of diffuser (or volute) blades, yielding a 
positive normalized precession frequency. 

f* = n x /n = /8|/(#V| - n 2 ) (27) 


Table 3 Dominant normalized precession frequency / + and peak-ex 
citation-amplitude frequency / * versus n, and n 2 ; iW 0) = 0.5 



na 

n 

r 

r 

(<7,»0.38 mm ) 

r 

( C p >0. T 3 mm ) 

8 

1 

7 

1.14 

4.2 

4.4 


2 

6 

1.33 

3.6 

3.8 


3 

5 

1.60 

2.9 

3.1 


4 

4 


2.1 

2.4 


5 

3 


1.4 

1.6 


6 

2 

4 

0.8 

0.9 


7 

1 

S 

0.3 

0.3 

7 

1 

6 

1.16 

3.6 

3.8 


2 

5 

1.40 

2.9 

3.1 


3 

4 

1.75 

2.1 

2.4 


4 

mm 

2.33 

1.4 

1.6 


5 

2 

3.5 

0.8 

0.9 


6 

1 

7.0 

0.3 

0.3 

6 

1 

5 

1.2 

2.9 

3.1 


2 

4 

1.5 

2.1 

2.4 


3 

3 

2.0 

1.4 

1.6 


4 

2 

3.0 

0.8 

0.9 


5 

i 

6.0 

0.3 

0.3 

5 

! i 

4 

1.25 

2.1 

2.4 


2 

3 

1.67 

1.4 

1.6 


3 

2 

2.5 

0.8 

0.9 


4 

1 

5.0 

0.3 

0.3 
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Fig. 9 I £,(0.95) I versus /; n * 3, u^O) = 0.5; C, = 0.36 and 0.72mm 


Moreover, for significant amplification within the leakage an- 
nulus, / + must lie near/*, the peak-amplitude excitation fre- 
quency. Table 3 shows the variation of/ + and /* for various 
combinations of n x and n 2 . The case of n x = 8 (eight-biaded 
impeller) and n 2 - 4 (four-bladed diffuser) yields a close prox- 
imity of /* = 2.1, 2.4 to / + = 2.0; however, this is an 
unrealistic combination. For a practical configuration, the 
nondimensional frequencies/ * = 1.67 and/* = 1. 6 are closest 
for ri\ = 5(five-bladed impeller) and n 2 = 2(double-discharge 
volute). Figure 9 illustrates I^K.95)! for 0) = 0.5 and n 
= 3, confirming the predictions of Table 3. An amplification 
by a factor of 2.6 is predicted for new clearances and 3.2 for 
worn clearances. 

Numerical Uncertainty. The numerical uncertainty issue 
for the results presented concerns the numerical integration of 
Eq. (17). The results presented were obtained with a fourth- 
order Runge-Kutta integrator package using 200 integration 
steps for the interval [0, 1]. Repeating these calculations with 
400 integration steps yielded the same results to about three 
significant figures. 
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Summary and Conclusions 

An analysis has been developed and results presented for 
the pressure oscillations in the leakage annulus between a 
shrouded pump impeller and its housing. These pressure os- 
cillations are driven by a circumferential variation of the im- 
peller discharge pressure which can precess either in the same 
or opposite to the direction of rotation. The circumferential 
variation can be modeled with a Fourier decomposition with 
each mode having n diametral nodes across the impeller. The 
peak-pressure oscillations within the impeller are predicted to 
occur near the exit wearing-ring seal in association with a 
centrifugal-acceleration-mode response of the fluid within the 
annulus (Childs, 1988, 1990b). The peak-amplitude-excitation 
frequency increases linearly with n. Using Bolleter’s (1988) 
work which provides a dependency of n and the precessional 
frequency on the number of impeller (n i) and diffuser ( n 2 ) 
blades, situations are presented which can yield large ampli- 
fications (or significant attenuation) of impeller discharge var- 
iations. The occurrence and nature of the pressure oscillations 
are shown to depend on: (a) the tangential-velocity ratio of 
the fluid entering the seal, (b) the order of the Fourier coef- 
ficient, (c) the closeness of the precessional frequency of the 
rotating pressure field to the first natural frequency of the fluid 
annulus, and (d) the clearance of the wearing-ring seal. 


The present results suggest an explanation for the internal 
melting observed on SSME HPFTP seal parts. However, given 
liquid hydrogen’s significant compressibility, a more complete 
analysis, including fluid compressibility, is in order. 
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Compressibility Effects on Rotor Forces in the Leakage Path between a Shrouded 
Pump Impeller and Its Housing. (August 1993) 

Nhai The Cao, B.S., Texas A&M University 


Chair of Advisory Committee; Dr. Dara Childs 


A modified approach to Childs' previous work (1989,1992) on fluid-structure 
interaction forces in the leakage path between an impeller shroud and its housing is 
presented in this paper. Three governing equations consisting of continuity, path- 
momentum, and circumferential-momentum equations were developed to describe the 
leakage path inside a pump impeller. Radial displacement perturbations were used to 
solve for radial and circumferential force coefficients. In addition, impeller-discharge 
pressure disturbances were used to obtain pressure oscillation responses due to 
precessing impeller pressure wave pattern. Childs' model was modified from an 
incompressible model to a compressible barotropic-fluid model (the density of the 
working fluid is a function of the pressure and a constant temperature only). Results 
obtained from this model yielded interaction forces for radial and circumferential force 
coefficients. Radial and circumferential forces define reaction forces within the 
impeller leakage path. 

An acoustic model for the same leakage path was also developed. The 
convective, Coriolis, and centrifugal acceleration terms are removed from the 
compressible model to obtain the acoustics model. The compressible model is 
compared with the incompressible model and the acoustic model. A solution due to 
impeller discharge pressure disturbances model was also developed for the 
compressible and acoustic models. The results from these modifications are used to 
determine what effects additional perturbation terms in the compressible model have 
on the acoustic model. 

The results show that the additional fluid mechanics terms in the compressible 



model cause resonances (peaks) in the force coefficient response curves. However, 
these peaks only occurred at high values of inlet circumferential velocity ratios, 
w eo(0)>0.7. The peak pressure oscillation was shown to occur at the wearing ring 
seal. Introduction of impeller discharge disturbances with n = 11 diametral nodes 
showed that maximum peak pressure oscillations occurred at nondimensional precession 
frequencies (J — Q/co where © is the running speed of the pump) of / = 6.4 and / 
= 7.8 for this particular pump. Bolleter's results suggest that for peak pressure 
oscillations to occur at the wearing ring seal, the nondimensional excitation frequency 
should be on the order of / = 2.182 for n = 11. The resonances found in this 
research do not match the excitation frequencies predicted by Bolleter. At the 
predicted peak excitation frequencies given by Bolleter, the compressible model shows 
an attenuation of the pressure oscillations at the seal exit. 

The compressibility of the fluid does not have a significant influence on the 
model at low values of nondimensional excitation frequency. At high values of 
nondimensional frequency, the effects of compressibility become more significant. For 
the acoustic analysis, the convective, Coriolis, and centrifugal acceleration terms do 
affect the results to a limited extent for precession excitation and to a large extent for 
a pressure excitation when the fluid operates at relatively high Mach numbers. 


To the memory of my grandfather, Cao Van Tac. 
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CHAPTER I 
INTRODUCTION 

In the past, wear-ring seals used on the SSME HPFTP (Space Shuttle Main 
Engine High Pressure Fuel Turbopump) made from KEL-F plastic came back after 
operation revealing highly unusual characteristics. Despite being in constant contact 
with liquid hydrogen, post-test inspection showed that interior points in the stator 
element had melted and resolidified. The material used in the seal stator has poor heat 
conduction properties and high internal hysteretic damping. Pressure oscillations 
adjacent to the seal may be a source of cyclic stress producing hysteretic losses. This 
investigation will examine the pressure oscillations which may cause cyclic stresses in 
the leakage path between the impeller shroud and its housing. Possible sources of 
excitation causing the seal to melt and resolidify will be investigated. 

For this project, an analysis will be performed for a bulk flow model of the 
leakage path between a pump impeller shroud and a housing along the front side of the 
impeller, from inlet to discharge (Fig. 1). Simply defined, a bulk flow model 
considers only the average of the velocity distribution across the flow field. The 
research will be an extension of analyses performed previously by Childs (1989, 1992) 
for a shrouded pump impeller and its housing. 

The working fluid will be modelled as a barotropic fluid in this analysis, instead 
of an incompressible fluid, to account for fluid compressibility. The density and 
viscosity of barotropic fluids depend only on the local pressure and are independent of 
temperature. This assumption is reasonable for most cryogenic fluids, where viscosity 
is low and effects of viscous heating are negligible. The properties of the working 
fluid, i.e., density and viscosity, will be implemented into a new analysis by using a 
general 32-term, thermodynamic, equations-of-state program, MIPROPS (McCarty, 
1986, modified by San Andres, 1991). 


This paper is modelled after the ASME Journal of Tribology. 



Figure 1 - Impeller stage and surface geometry 
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CHAPTER n 
LITERATURE REVIEW 

Childs (1989) performed a bulk-flow analysis for the leakage path between an 
impeller shroud and a pump housing. Three governing equations of motion were 
derived for a bulk-flow model to represent incompressible fluid flow in the leakage 
path of a conventional water pump impeller. Three equations, consisting of a 
continuity equation, a path-momentum equation, and a circumferential-momentum 
equation, were used to solve for rotordynamic forces due to a precessional excitation 
of the rotor. 

Childs used a perturbation expansion in the eccentricity ratio of the governing 
equations of fluid motion for small motions about a centered impeller position yielding 
a set of zeroth and first-order governing equations. A zeroth-order solution was 
obtained by an iterative procedure to define the leakage, pressure, and circumferential- 
velocity distribution. Using a perturbed clearance function due to a radial displacement 
perturbation, Childs evaluated the first-order model at several inlet circumferential 
velocity conditions to obtain the first-order perturbed solutions. First-order 
perturbation results provided rotordynamic coefficients (direct and cross-coupled 
stiffness, damping, and mass) and lateral reaction forces for the model. Childs' 
predictions for the impeller of Fig. 2 are shown in Fig. 3. The predicted radial and 
circumferential force coefficients are shown versus the nondimensional precessional 
frequency for nondimensionalized inlet circumferential velocities of u eo (0) = 0.5, 0.6, 
and 0.7. Nondimensional precessional frequency is the ratio of the impeller precession 
frequency, 12, to its running speed, to. The graphs showed a considerable dip , or 
resonance, in the radial and circumferential force response coefficients at higher values 
of Ue 0 (O). The radial and circumferential force coefficients represent the 
nondimensionalized reaction forces acting on the impeller face due to impeller 
precession. Childs showed that the centrifugal acceleration terms in the momentum 
equations produced the "dip" in the results. By removing the 
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DIMENSIONS in 
62 — MILLIMETERS 



Figure 2 - Nominal configuration of conventional water impeller 




Figure 3 - Nondimensional force coefficients for the conventional impeller (a) 
circumferential-force coefficient, (b) radial-force coefficient 
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centrifugal acceleration term from the path-momentum equation, the "dips" in the plots 
were eliminated. 

Bolleter (1988) presents a relationship between the difference ( n ) in the number 
of impeller blades (n,) and diffuser vanes (n 2 ) and the precessional frequency for 
pressure pulsations in an impeller leakage path. Various combinations of impeller 
blade and diffuser vane number causing vibrations and pressure pulsations in the 
impeller are described. The relationship presented by Bolleter states that a pressure 
pattern develops with n — I - n 2 I diametral nodes around the impeller exit. The 
precession velocity of the pressure oscillation is /!©//( I n t - n 2 b- F° r example, an 
impeller with n t — 1 1 blades and n 2 — 6 diffuser vanes would have amplifications of 
the pressure oscillations at frequencies of multiples of ll©/(lll - 6l)=2.2co. 

Childs (1992) performed an analysis similar to (Childs, 1989) incorporating the 
effects of different numbers of impeller blades and diffuser vanes for the bulk flow 
model using excitations due to discharge-pressure oscillations instead of orbital motion. 
This analysis considered the harmonic response of flow within the annulus due to 
variations in the discharge pressure of the impeller. This analysis also compared the 
effect on the response of the pressure oscillation due to different numbers of pump 
impeller blades and diffuser vanes. Zeroth and first-order perturbation equations were 
also derived for this analysis. However, the first order perturbations in this analysis 
were excited by discharge-pressure perturbations instead of impeller precession. The 
impeller discharge excitation was defined as a precessing harmonic pressure oscillation 
with n nodes and a precessional frequency of Ci. 

Results from the analysis due to perturbed flow in the leakage path caused by 
oscillations in the impeller discharge pressure show that the peak pressure oscillation 
occur near the exit ring seal. The pressure oscillations from the impeller were shown 
to depend on the circumferential velocity of the fluid entering the seal, the Fourier 
coefficient, n, and the relative closeness to the first resonant frequency of the fluid to 
the peak precessional frequency of the rotating pressure field. Note that n represents 
both the Fourier coefficient used by Childs and Bolleter's n = \n t - n 2 \. 


6 


Acoustic modes are produced by the interaction of fluid inertia and 
compressibility. Thompson (1988) explains that the wave equation, the fundamental 
equation of acoustics, is obtained by assuming that the convective acceleration terms 
are negligible compared with the temporal acceleration terms. For ordinary acoustic 
analysis, Thompson states that this assumption is "highly satisfactory" for fluid flow 
characterized by a low Mach number, typically M 2 « 0. 1 . By removing the convective, 
Coriolis, and centrifugal acceleration terms from a modified general perturbation 
(compressible) version of Childs' model, a similar wave equation can be obtained. 
This resulting equation, in theory, can be used for an "acoustic" analysis of the flow 
fields. 

San Andres (1991) developed a solution procedure for a model of fluid flow in 
turbulent hydrostatic bearings and annular seals operating with cryogenic barotropic 
fluids. He used a 32-term equations-of-state program provided by NBS Standard 
Reference Data Base for prediction of the properties of LH 2 , L0 2 , LN 2 , and other 
fluids at different pressures and temperatures. The code, MIPROPS, delivers fluid 
properties which are used in the analysis procedure. In addition to obtaining the fluid 
properties from MIPROPS for use in a compressible model, San Andres also 
considered the properties of the working fluid as a linear function of pressure. From 
his results, San Andres found that for highly compressible fluids, such as liquid 
hydrogen, the barotropic properties model based on an equation of state gave accurate 
leakage and force response for bearings and seals with a large pressure differential. 
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CHAPTER m 
OBJECTIVES 

This research will introduce compressibility of the working fluid in the leakage 
path as an extension to the analyses performed by Childs (1989, 1992). The results 
will provide information concerning the relationship between an incompressible, a 
compressible, and an acoustic model with liquid hydrogen as the working fluid. The 
results will also provide information regarding the cause of the unusual behavior 
exhibited by the KEL-F plastic rotor element and verify the validity of Thompson’s 
assessment concerning the effects of fluid mechanics and acoustics. This research 
project will also analyze the effects of the centrifugal acceleration modes and acoustic 
modes of a barotropic fluid in the leakage path between a shrouded pump impeller and 
its housing. 

The results obtained from a compressible-flow model will be compared with 
the results of an incompressible model. The compressible-flow model will also be 
reduced to an acoustics model, the results of which will be compared to the 
compressible-flow model. This comparison will be performed for two different 
perturbation excitations: (a) a precessional excitation involving an orbital motion of 
the rotor, and (b) a pressure oscillation excitation, involving perturbation of the 
discharge or inlet pressure of the leakage path. The geometric and operating 
characteristics of the first stage impeller of the SSME HPFTP will provide the 
parameters used for the governing equations. 

A bulk-flow model will be developed and used to simulate the leakage path 
inside the first impeller stage of the Space Shuttle Main Engine Turbopump. Results 
from the computation should indicate if any interaction exists between acoustic and 
centrifugal acceleration modes, and the influence of fluid mechanics terms (convective, 
Coriolis, or centrifugal acceleration terms) on acoustic modes. 
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CHAPTER IV 

GEOMETRIC AND OPERATING CHARACTERISTICS 

Figure 4 shows the first-stage impeller of the Space Shuttle Main Engine 
Turbopump. The impeller measures 0.3048 m in diameter at the exit (leakage path 
inlet) and 0. 1905 m in diameter at the seal inlet (leakage path exit). The impeller is 
also characterized by 24 impeller blades and 13 diffuser vanes. At full power levels, 
this stage operates at 34,000 rpm, with an inlet pressure at the entrance to the leakage 
path of 13.79 MPa and a discharge pressure of 1.72 MPa at the exit of the wearing 
ring seal. Operating at a pump speed of 34,000 rpm, the resultant velocity vector of 
the liquid hydrogen inside the leakage path is calculated to be about 0.4 times the 
acoustic velocity of liquid hydrogen. 

The wearing ring seal contains four steps, measuring 0.1915 , 0.185 , 0.1786, 
and 0.172 m in diameter, which accommodate four teeth at the end of the impeller 
blade. The radial clearances between the seal and the rotor teeth are estimated to be 
0.229 mm (0.009 in), accounting for radial expansion of the rotor during operation. 
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Figure 4 - SSME HPFTP first impeller stage 
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CHAPTER V 

BULK FLOW MATHEMATICAL MODEL 


Childs’ governing equations will be modified for this project to reflect a bulk 
flow model operating with a barotropic fluid. As in Childs’ analysis, these equations 
will be nondimensionalized and perturbed to yield zeroth and first-order governing 
equations. 


5.1 General Governing Equations 

Using the approach taken by Childs (1989), the governing equations are: 
• Continuity Equation 


dpH 

dt 




. "»pI7. - 0 

RdS s 


( 1 ) 


Path-Momentum Equation 
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• Circumferential-Momentum Equation 
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R ae 
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An additional governing equation is obtained by using MIPROPS to obtain the 
properties of liquid hydrogen. 



• Equations of State 
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p=p ( pj °) 

P=p(/>,r“) 


(4) 


With the exception of the addition of p to the continuity equation, Eqs. (1-3) 
are identical to Childs' (1989). The path and circumferential-momentum equations do 
not change from Childs' model because the density (p) drops out of the momentum 
equations when the continuity equation is used to simplify them. 

The equations of state in Eq. (4) define the density and viscosity for the bulk- 
flow model. The variation in density of the fluid in the model will be implemented by 
assuming that the working fluid is barotropic. Here, the variation in density of the 
working fluid will be modelled as a function of pressure and a constant temperature 
only. In this investigation, the MIPROPS code calculates the value of density and 
viscosity at a constant temperature of 23.37 K with varying input pressures, and 
returns values of density and viscosity to the main program. The variation of viscosity 
with respect to pressure were very slight; therefore, viscosity was kept constant. 

H(5, 0, t) in the governing equations defines the clearance between the impeller 
and the housing. Nondimensionalization of this variable is given in the nomenclature 
and also later in this text. 

Hir's (1973) definitions were used to define the shear stress components of the 
rotor and stator surfaces. The equations shown below define the shear stress acting 
on the impeller and its housing. The first subscript in the equations denotes the 
direction of fluid flow (path and circumferential), and the second subscript refers to the 
surface (stator and rotor), respectively. 


JHJ+1 

t. ■ fp<X"[i *( u J u ,f \ 2 


( 5 ) 


12 


mr* 1 

- f PWa ^ 1 + (W 2 ] 2 


( 6 ) 

(7) 


mr+1 

v - ^i‘W>-* a K{ i *w>- Ra y v fi 2 


( 8 ) 


Reynolds' number used in these equations is represented by, 

R a = 2HUJv 


(9) 


Boundary Conditions 

The pressure drop of the inlet to the leakage path provides the inlet boundary 
condition (s=0) given by the relationship, 

p J -p 0 (o,e,i)-p(i»5)(/>,e,!)/2 <1#) 

The exit wearing-ring seal defines the following exit boundary condition, 




(ID 


These boundary conditions apply directly for precession excitations. 

For the analysis which examines the changes in (u } , u e , p) due to changes in 
the impeller's discharge P t or inlet pressure P s , the following boundary conditions are 
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stated for the inlet and exit, respectively, 

PjfiA - Pfi M * p(i (12) 

P(t,,e,() - p,(0, r) - 1 c+ufaej) (13) 

Eqs. (12-13) differ from Eqs. (10-11) because P e and P s are now also functions of 
time. 


5.2 General Perturbed Equations - Nondimensionalization and Perturbation 
Analysis 

Introducing the following variables into Eqs. (1-8), 

u,-UJV, « e = UJRp, p = P/pVf, p=p/ Pi 

h=H/C p s = S/L s , r = m i (14) 

t=gU, i> = F/ZfjCo, T=LJV l 


yields nondimensional governing equations. 

The perturbation variables used to obtain zeroth and first-order equations are 
defined by, 

u s~ u so + eu si> P = Po +€ Pi (15) 

«e = «eo + €W ei> P = Po + *Pi 

where e is the perturbation coefficient to be defined separately below for precession 
and discharge-pressure excitation. 

5.2.1 Zeroth Order Solution 

The zeroth-order equations are the same for the precession and exit-pressure 


excitations. 


14 


5.2. 1.1 Zeroth-Order Equations 

The path and circumferential velocity distribution and the leakage rate for a 
centered impeller position are defined by the following zeroth-order governing 
equations. 


• Continuity Equation 


rh o u *oPo " 1 


• Path-Momentum Equation 


l_dPo 

Po 


du 


=u 


so 


sO 


ds 


1 dr 
r ds 



+ 


K + °r) 
2 



• Circumferential-Momentum Equation 
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• Equation of State 

P = p(pj°) 


( 16 ) 


( 17 ) 


( 18 ) 


( 19 ) 


( 20 ) 


Eq. (19) represents friction factor definitions for the stator and rotor surfaces, 
respectively, introduced in Eqs. (5-8). 
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5.2. 1.2 Zeroth-Order Solution 

The zeroth-order continuity Eq. (16) can also be expressed as, 

1 1 dpp + l 

h 0 ds p 0 ds r ds 


du 




ds 


= -u 


sO 


This equation can be substituted into Eq. (17) to obtain 


1 dp 0 2 
— — 1/ 

1 dh^ 1 dp 0 l dr 

l \ 2 

_ 1 drl w eo 

* » ■ sO 

Po 

h 0 ds p 0 ds r ds 

r ds\ b j 


(°s +a r )„2 
^ U sO 


( 22 ) 


The governing zeroth-order equations now reduce to two governing equations, 
consisting of Eq. (18) and Eq. (22). Eq. (20) defines p solely as a function of p and 
provides the density used in Eq. (18) and Eq. (22). 


Boundary Conditions 

The inlet boundary condition for the zeroth-order pressure relationship can be 
expressed from Eq. (10) as, 


P Q (0) = 




(23) 


The zeroth-order solutions are obtained by solving Equations (18) and (22) 
iteratively. An initial (s=0) fluid velocity V; is estimated which then defines ujs). 
A specified u M (0) and the calculated p 0 from Eq. (23) are used to numerically integrate 
the zeroth-order equations (18) and (22) from the path entrance (s=0), to the path exit 
(s=l). The procedure is continued with revised values of V ( until convergence is 
obtained between the prescribed and the calculated exit pressure. 


5.2.1. 3 Zeroth-Order Results 

The zeroth-order results provide the mass flow rate through the leakage path. 
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The discharge coefficient C dt used in this analysis was obtained by using the leakage 
rate through the wearing ring seal. The leakage through the wearing ring seal of the 
SSME HPFTP was calculated using a seal leakage code developed by Morrison et al. 
(1983), and this value was used to calculate the discharge coefficient at specific 
impeller operating conditions. The seal leakage code uses the geometry of the seal 
and the operating conditions, i.e., inlet and exit pressures, viscosity, density, etc. to 
calculate the leakage rate through the seal. Note that this seal leakage code treats the 
working fluid as incompressible. 

The seal leakage code yielded several mass flow rates for the prescribed 
geometry and operating conditions, depending upon the pressure drop across the 
impeller and the wearing ring seal. Several discharge coefficients were tested in the 
model to match the flow rate through the impeller leakage path and the flow rate 
through the seal. Once the two flow rates converged, the resulting C& was used in 
the model as an exit restriction boundary condition. The mass flow rate through the 
seal was found to be 1.6373 kg/s, with the resulting C de being 7039.6, and inlet and 
exit pressure of the seal being 8.101 MPa and 1.72 MPa, respectively. This pressure 
drop across the seal represents about one third of the total pressure drop across the 
entire impeller leakage path. 

Zeroth -order pressure solutions for the incompressible and compressible models 
are shown in Figure 5. The nondimensional path velocity along the leakage path is 
given in Figure 6, and the zeroth-order circumferential path velocity is shown in 
Figure 7. For the incompressible model, a mean value for the density of liquid 
hydrogen inside the impeller leakage path was used to obtain the results shown. 

The results of the compressible model vary only slightly from the 
incompressible model for the zeroth-order solution. The pressure distribution across 
the leakage path shows the same trend and approximately the same magnitude of 
pressure drop, but not exactly the same inlet and exit pressure values. The C* found 
earlier provides the exit restriction for the seal and therefore is used for both models. 
Because the incompressible model uses an average density along the leakage path, the 
inlet and exit densities for the two models will be slightly different. This accounts for 



Nondimensional Path Length (s) 


Figure 5 - Zeroth-order pressure p distribution along leakage path for compressible 
and incompressible models for u^O) = 0.7 



Nondimensional Path Velocity (us(0)) 
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Figure 6 - Zeroth-order path velocity (u^) distribution along leakage path for 
compressible and incompressible models for u^O) = 0.7 


mensionai Circumferential velocity lutn^u^ 
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the difference in the dimensionless pressure magnitudes shown in Figure 5. The 
dimensional pressures at the inlet are almost exactly the same value, confirming the 
inlet boundary condition. As shown in Figure 6 and Figure 7, the zeroth-order path 
velocity and the zeroth-order circumferential velocity are not affected significantly by 
compressibility. 


5.2.2 First-Order Equations 

First-order governing equations define the path and circumferential velocity and 
pressure distribution along the leakage path due to perturbed clearance function or 
perturbed discharge-pressure. 

First-order perturbation equations obtained by the perturbation expansion of Eq. 
(14) are 

• Continuity Equation 
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• Circumferential-Momentum Equation 
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• Equation of State 

dp, dp l dp l dp 1 _dp ± dp 1 dp, s dp 1 ty 1 (27) 

1 T = ’ dx = d P dx’ ae rfp so 

With the exceptions of A 4> and A 4 „ which are defined in the appendix, the parameters 
A,„ A 2j , etc., in these equations can be found in Childs (1989). 

The dependency of p with respect to 5, r, and 6 in Eq. (24-26) are eliminated 
from the governing equations by applying the definitions of Eq. (27). The relationship 
between p and p in Eq. (27) was obtained from the results produced by MIPROPS. 



22 


CHAPTER VI 

FIRST ORDER EQUATIONS AND SOLUTIONS FOR GENERAL 
PERTURBED COMPRESSIBLE MODELS 

This chapter provides the first order equations and solutions for precession 
excitation and pressure excitation cases for the general perturbation (compressible) 
models. Results for both excitation cases are presented using an inlet tangential 
velocity of u eo (0) = 0.7. The radial and circumferential-force coefficient response 
curves represent the nondimensional reaction force acting on the impeller face in the 
respective directions versus nondimensional frequency ratio, /. 

For the first-order solution analysis, a separation of variable approach was used 
to obtain complex ordinary differential equations. The resulting coupled equations 
were integrated to obtain the nondimensional radial and circumferential force response 
coefficients. The calculated results from the first-order precession excitation and the 
pressure excitation of the general perturbed bulk flow model provide predictions which 
can be used to qualify and quantify the effects of fluid compressibility in the model. 

6.1 Precession Excitation: General-Perturbation Model 

The precession excitation of the general perturbed model uses the perturbed 
clearance function, 

= A, c &T)cose + /t u (s,r)sin0 ( 28 ) 

as the excitation. 

The theta dependency of Eqs. (24-26) can be eliminated by substituting the 
following solution format 

u sl = K, ic cos0 + « jif sin0 « e , = w 01f cos0 + « 6Lt sin0 

Pi = Pie cos 0 + /> Ls sin0 p , = £ lc cos0 + p ^sin© 

into Eqs. (24-26), which yields six real equations. 
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Three complex equations in the independent variables s and x can be obtained 
by introducing the complex variables 

U s r u sic + J u sis> =M eic + 7W = Pic + ./Pu (30) 

=Pic+JPw A 1 “* le + A* 


These complex equations in the independent variables s and x are 
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Since the equation of state is a function of pressure and a constant temperature only, 
it can be modified as shown in Eq. (27) and used to remove the dependency of p, from 
s in Eq. (31). Further simplification of Eqs.(33-35) can be made by using the 
following definitions provided by Childs (1989) 



(34) 



( 35 ) 
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where 


q=x + jy 


(36) 


Eq. (36) represents the physical motion of the rotor in the x and y directions, shown 
in Figure 1. 

Assuming a harmonic seal motion of the form 

<7=g 0 ^’ / =Q /<*> (37) 

the corresponding harmonic solutions can be stated, 

* s} " V*. «e, =^ 1 ^ (38) 

which yields the following three complex ordinary differential equations of motion, 
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Elements used in Eqs. (40-41) can be found in Appendix A. 
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6.1.1 Boundary Conditions 

Nondimensional first-order boundary conditions for the precession excitation can 
be stated from equations (10) and (11) as 
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Additionally, the perturbation entrance circumferential velocity can be stated as zero, 

5,,(0)=0 (44) 

The solution to this set of equations is obtained by applying the procedure 
presented by Childs (1989). The solution procedure used to determine the reaction 
forces and moments are also given by Childs (1989). 


6.1.2 First-Order Results 

Radial and circumferential force response coefficients for general compressible 
and incompressible precession excitation models operating with an inlet circumferential 
velocity of u^O) = 0.7 are shown in Figure 8 and Figure 9. The results 



Nondimensional Frequency (f) 

Figure 8 - Radial response for compressible and incompressible models for precession 
excitation for uU 0) =0.7 






Nondimensional Force Response 
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Nondimensional Frequency (f) 


Figure 9 - Circumferential response for compressible and incompressible models for 
precession excitation for n 0o (O) = 0.7 
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show that the effects of compressibility do not affect the results of the model to a 
great extent over the frequency range considered. As with Childs’ results, local 
resonance peaks occur in the response curves. Two local peaks occur in the response 
curves, at nondimensional frequency ratios of / = 0.1 and / = 0.8. For low inlet 
values of circumferential velocity w^O) = 0.5, the behavior of the response curves 
of the compressible model is virtually the same as the results shown for an 
incompressible model shown by Childs (1989), i.e., the resonance in the response 
curves diminish at lower values of inlet circumferential velocity. 

Complex first-order pressure, path velocity, and circumferential velocity for 
different nondimensional excitation frequencies along the leakage path produced from 
the first-order perturbation analysis provide approximate complex modes at the 
resonant frequencies. The real and imaginary parts of these results are used to obtain 
amplitude and phase plots at the resonant frequencies. Amplitude and phase plots of 
first-order nondimensional pressure in the leakage path at the frequency ratios of / = 
0.1 and / = 0.8, where the local peaks occur, are shown in Figures 10 and 11. 
Complex modes for u 5l at the same frequency ratios are shown in Figures 12 and 13. 

At the frequency ratio of / = 0.1, the amplitude of the pressure along the 
leakage path steadily increases and has a maximum value at the exit of the leakage 
path. For the frequency ratio of / = 0.8, the maximum pressure amplitude occurs 
near the middle of the leakage path. The mode shape for u sl at / = 0.1 shows a 
slightly decreasing amplitude along the leakage path, with a minor increase in the 
amplitude at the exit. As with the mode shape for the pressure distribution at / = 
0.8, the maximum magnitude occurs near the middle of the seal. 

As with Childs’ model, the resonant peaks found in this analysis can be 
attributed to the centrifugal acceleration terms. When the centrifugal acceleration 
terms are removed from the model, the local peaks in the response curves diminish 
significantly. The mode shapes at / = 0.8 do not support the theory that maximum 
pressure oscillations at the leakage path exit are causing the seal to melt because the 
maximum amplification occurs near the middle of the leakage path. 
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Figure 10 - Amplitude and phase plot of nondimensional pressure p for precessional 
excitation of compressible model at / = 0.1 
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Figure li - Amplitude and phase plot of nondimensional pressure p for precessional 
excitation of compressible model at / = 0.8 
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Figure 12 - Amplitude and phase plot for precessional excitation for path velocity u s2 
of compressible model at / = 0.1 
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Noticeable effects of compressibility upon the model can be seen at higher 
nondimensional frequency ratios. Figure 14 shows the response curves of the 
compressible and incompressible models for high frequency ratios. The force 
responses of the incompressible model continue to grow with increasing/, while the 
results of the compressible model exhibit peaks corresponding to the natural 
frequency. This is consistent with classical vibration analysis, with the incompressible 
model having no natural frequency due to an infinitely high bulk modulus, while the 
compressible model has a natural frequency corresponding to the compressibility of 
liquid hydrogen in the model. 

6.2 Pressure Excitation 

This part of the research involves introducing a time and circumferentially 
varying impeller discharge pressure to the compressible model. With different 
excitations involved in this model, new boundary conditions also exist. The pressure 
perturbation takes the form of the following equation. 

• Supply Pressure-Excitation Perturbation 

(4S) 

The occurrence of epsilon in this equation represents the perturbation coefficient for 
discharge-pressure excitation. 

6.2.1 First-Order Equations 

Nondimensionalization and perturbation of the general governing equations (6- 
8) yields the same governing equations for pressure excitation perturbation as for 
precession excitation governing Eqs. (31-33), except with h, — 0 for this model. 
Therefore, the first-order governing equations for pressure excitation will not be 
presented in this section. 



Nondimensional Frequency (f) 


Figure 14 - Nondimensional force response of compressible and incompressible models 
for precessional excitation for extended frequency range 
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6.2.2 Boundary Conditions 

Nondimensionalization and perturbation of the boundary conditions introduced 
in Eq. (12-13), the first-order boundary condition can be stated as 

p,m -p,mt) - (i *5) U „(o,e,fl|i io) 

pflfiA -p.m * W !>) 

p sl (d,t), Ptift’O provide the boundary excitations and take the form, 

PsiM * e jDl (p slc cosnd +p sl ?mnQ) (47) 

p ei m=e ja, lj)' ]c cosnQ *p tl ^nnQ) 


6.2.3 First Order Solution 

The theta and time dependency of the first-order pressure excitation governing 
equations is eliminated by assuming 

u s i =e^(u sU cosnQ + u sI ^mnQ) 

« e j * (u 61c cos nd + « fiJ sinn0) 

p x =e jft (p lc cosnQ +/> u sinn0) 


where n represents the difference of number of impeller blades and vanes, also defined 
as the number of diametral nodes (Bolleter). Substitution of Eqs. (48) into the first- 
order pressure excitation governing equations, equating coefficients of cos nd and sin 
nd, and using complex variables described in Eq. (30) reduces the real equations to 
three, complex ordinary differential equations 


d_ 

ds 


U sl 



«ei 

>+[B(nJ,s)) 

«ei 

.?«] 


}\ 


( 49 ) 
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where 
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•^29 +JTT 

A 

dp 0 . bL 

— i 

(50) 


U s0 U s0 

46 ds J r Us0 p 0 *. 


*31 

B 32 


-^4 s _ *33 dPo 

P 0 P 0 * 



Eq. (49) differs from Eq. (40) in that Eq. (49) is now a function of n with h, = 0. 

These equations are solved using the method presented by Childs (1992). 

6.2.4 First-Order Results 

Results from the pressure excitation perturbation analysis using a difference of 
impeller blades and vanes (n) is presented in this section. Since the SSME HPFTP 
wearing-ring seal is located at the end of the leakage path (s = 1), the results from 
this section will concentrate on this area of the leakage path. The pressure-oscillation 
response at the end of the leakage path versus nondimensional excitation frequency 
corresponding for n = 11 is illustrated in Figure 15. The value of n = 11 is used 
because, as mentioned previously, the SSME HPFTP has 24 impeller blades and 13 
diffuser vanes, for a difference of 11. Bolleter (1988) explains that the dominant 
pressure pattern to be expected for this arrangement would have 1 1 diametral nodes. 
This is important because it is the value used as n in Eq. (48). Figure 15 shows that 
the peak pressure oscillation at the exit of the impeller leakage path (5=1) occurs at 
a nondimensional frequency of about / = 6.4 and / = 7.8. Bolleter predicts that the 
peak pressure oscillation should occur as a function of the difference in the number 
of impeller blades and vanes, presented in Tables 1 and 2. 



38 


The resonant frequency ratios at / = 6.5 and / = 7.8 clearly are not functions 
near the predicted excitation frequencies provided in Table 2 for n = 11. None of the 
other cases tested yielded results which correspond to the values shown in Table 2 for 
the respective value of n. These results reveal that amplification of peak pressure 
oscillation does not occur at the expected frequency ratio. At the predicted frequency 
ratios, the cases tested for this model showed that an attenuation, rather than an 
amplification of the pressure waves occurred. 


Table 1 - (n, - for various combinations of multiples of impeller blades and vanes 


n, 

24 

48 

96 

13 

11 

35 

83 

26 

-2 

22 

70 

39 

-15 

9 

57 


Table 2 - Expected peak nondimensional frequency ratios, / or, {n/(n, - n^}, for 
various combinations of multiples of impeller blades and vanes 


n 2 n, 

24 

48 

96 

13 

2.18 

1.37 

1.16 

26 

-12.00 

2.18 

1.37 

39 

-1.60 

5.33 

1.68 


Amplitude and phase plots (Figures 16 and 17) of the pressure oscillation for 
the two resonant frequency ratios found for n for the compressible pressure excitation 
model show that the peak pressure oscillation occurs at the exit of the leakage path. 
This result helps to support the hypothesis that pressure oscillations at the leakage path 
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Figure 16 - Amplitude and phase plot for pressure oscillation of compressible model 
inside impeller leakage path at / = 6.4 for n = 11 and u^O) = 0.7 






40 


Mode Shape 



Nondimensional Path Length (s) 


Phase Plot 



Nondimensional Path Length (s) 


Figure 17 - Amplitude and phase plot for pressure oscillation of compressible model 
inside impeller leakage path at / = 7.5 for n = 11 and u,( 0) = 0.7 




41 


exit are causing the seal to melt. However, the result does not provide enough 
evidence to provide a very plausible explanation for the behavior exhibited by the 
wearing ring seal. 

The response of the compressible model for other values of n is plotted in 
Figure 18. The graph shows the relationship of the peak pressure oscillation response 
to the number of diametral nodes. The nondimensional excitation frequency ratio at 
which the peak pressure oscillation occurs increases as the number of diametral nodes 
increase, consistent with Childs’ analysis (1992). However, unlike Childs results, 
where the magnitude of the peak amplitude stays relatively constant, regardless of n, 
the results produced from the compressible model show that the amplitude of the 
pressure oscillations increases as n increases. The effects of added compressibility to 
the results of this model is greatly enhanced compared to results of the precession 
excitation model. But, as with the precession excitation analysis, the influence due to 
compressibility effects is only noticeable at higher nondimensional frequencies. 
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Figure 18 - Pressure oscillation response for pressure oscillation of compressible 
model with multiple n and 0) = 0.7 


CHAPTER Vn 

FIRST ORDER EQUATIONS AND SOLUTIONS FOR ACOUSTIC MODELS 


For acoustics analysis, the temporal acceleration terms are the only acceleration 
terms which remain from the set of general governing equations. A highly reasonable 
assumption may be made that the contributions from the convective, Coriolis, and 
centrifugal acceleration terms are negligible compared to the temporal acceleration 
terms in the realm of ordinary acoustics (M 2 « 1) (Thompson, 1988). However, due 
to the high Mach number with which the working fluid in the SSME HPFTP operates 
(M = 0.4), the general perturbation fluid model and the acoustic model must be 
compared to investigate the effects of the convective, Coriolis, and centrifugal 
acceleration terms in an acoustic analysis. 

The zeroth-order solution for the acoustics models do not change from the 
previous results obtained because the zeroth-order equations and solution remain the 
same . 

7.1 Precession Excitation Model 

7.1.1 First Order Equations 

First-order governing equations for the acoustic precession excitation model can 
be obtained by removing the effects of fluid mechanics, i.e., the convective, the 
Coriolis, and the centrifugal acceleration terms from the general perturbation governing 
equations given in 5.2.2, Eqs. (24-26). Removing the convective acceleration terms, 
u^dug/ds ujiujds, the Coriolis acceleration term, lu^r/rds, and -2 u 0( flr/(b 2 rds), 
representing centrifugal acceleration, from the continuity Eq. (24) yields the first-order 
governing acoustic continuity governing Eq. (51). The first-order governing acoustic 
equations for the path and circumferential momentum are obtained by removing the 
convective acceleration terms from Eqs. (25-26). 
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Continuity Equation 


- dh i . a Pi A 1 
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dx 
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1 51s ds 1 ds 
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(51) 


Path-Momentum Equation 


1 dp, du,, . 

1 + w ei^2* + “A + Pi A « + T <*-Z- =h i A is 
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• Circumferential-Momentum Equation 

b h d P\ 


du 


~r~Rf> — 50^ + Wet ^ 20 + Wf/ ^ 3e + P^ 46 + 


(52) 


(53) 


7.1.2 First-Order Precession Excitation Solution 

The theta dependency of Eqs. (51-53) can be eliminated by substituting the 
solution format presented in Eq. (29). The three resulting complex equations in the 
independent variables s and r after introducing the complex variables of Eq. (32) are 
• Continuity Equation 


du 


Si _ juT 


ds 




1 1 dp c 


+ £j 
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.0 >r Wrfdpo 

J «An + — 

V“ Po 2& 


ds 


Po & 

K ^ 


= k 


u s0 dh 0 




.cor 

+ J-r~ u *> 

V 


“,0^1 cor^i ^i 

*o as p 0 4p, dr 


(54) 




Path-Momentum Equation 


d£. du 

*W< W =t,A 


U 


( 55 ) 


• Circumferential-Momentum Equation 


r R i Po 


+ a ei^e 



Vie 


(56) 


Assuming the harmonic seal motion described in section 6.1, Eq. (38), the three 
governing equations reduce to two differential equations of motion of the form, 


-M 

ds\p 


+[C] 




(57) 


where 



«£ o£Po r +F .IT<fe jUT ^i6 

Po dp q 2 J h 0 ds J r p 2i >+joTJ){ 

Q ^ 2^16 

^ (^2 


(58) 


The [C] matrix elements are given in Appendix C. 

The du ti /ds term has dropped from the governing equations and thus u eI can 
be solved directly without integrating the partial differential equations, yielding, 


‘ei 


Vie - 


±J*.-2!aX 

r Vo d P ) . 


l^n+juTf) 


(59) 


Results to these equations are obtained using the same solution procedure and 
boundary conditions described for the first-order precession excitation equations in 
section 6.2.2. The solution to these sets of equations can then be used to calculate the 
force and moment coefficients in the impeller leakage path for the acoustic case where 
classical acoustic assumptions are made. 
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7.1.3 First-Order Results 

The nondimensional radial force coefficient responses for the general 
perturbation and acoustic cases are shown in Figure 19. The response curves for the 
circumferential force coefficients are shown in Figure 20. The response curves for 
the acoustic case are less erratic in behavior, and the peaks exhibited by the 
compressible model totally disappear. This shows that the convective, Coriolis, and 
centrifugal acceleration terms do affect the results of the compressible model to a 
considerable extent at low frequencies. Childs (1989) correctly predicted that the 
centrifugal acceleration terms caused the resonance in the response curve. Although 
the results of the two models do not differ quantitatively, the trend exhibited by both 
models are similar, as evidenced by Figures 20 and 21. 

At higher values of nondimensional frequency, the effects of removing the fluid 
mechanics terms from the model is shown in Figure (21). The two models show 
about the same results. The natural frequency of the acoustic model is slightly lower 
and the peaks are slightly higher. Mode shapes for pressure oscillations at f — 11.5 
f — 12 are shown in Figures 22 and 23. The mode shapes show that, at these 
frequencies, the pressure oscillations are not the cause of the uncharacteristic behavior 
of the exit wearing ring seal due to a precession excitation. The mode shapes do not 
show conclusively that the pressure oscillations are occurring at the exit of the seal. 

7.2 Pressure Excitation Model 

7.2.1 First-Order Equations 

The equations obtained for the first-order acoustic pressure excitation 
perturbation use the same procedure as for the first-order acoustic precession 
excitation shown in section 7.1.1. The boundary conditions and perturbation 
excitation are defined by the equations used for the general first-order pressure 
excitation solution given in section 6.2. The resulting governing equations for 
pressure excitation of the acoustic model are the same as those given in section 7.1.1, 
except with h, being zero, analogous with the results obtained for the pressure 
excitation equations for the general perturbation model. 
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Nondimensional Frequency (f) 


Figure 19 - Radial response for acoustic and compressible models for precessional 
excitation for uJi 0) = 0.7 
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Figure 20 - Circumferential force response for acoustic and compressible models for 
processional excitation for uji 0) = 0.7 
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Figure 21 - Force response coefficient for acoustic and compressible models at 
extended frequency ratios 
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7.2.2 First-Order Pressure F.xcitation Solution 

Using Eq. (48) to eliminate the theta and time dependency from the first-order 
pressure excitation of the acoustic model and using the pressure excitation perturbation 
described in Eq. (45), the resulting governing equations can be represented as, 


Where elements of [D] are provided in Appendix D. 

In addition to h, being zero, Eq. (60) also is a function of n, unlike Eq. (57). 

7.2.3 first-Order Results 

Comparison of the general perturbation model and the acoustic model for n= 11 
shows that the peak pressure amplification around / = 7.8 displayed by the general 
perturbation model almost totally disappears for the acoustics model. From Figure 24, 
the results for the acoustics model show that no amplification of the pressure 
oscillation occurs, but rather, an attenuation of the pressure oscillation at all frequency 
ratios. It can be concluded from this result that the convective, Conolis, and 
centrifugal acceleration terms contribute significantly to the results of a bulk flow 
model when a pressure excitation is used as the perturbation parameter of the bulk 
flow model. The results here also differ with respect to the location of the resonant 
peaks in the nondimensional frequency range tested. For the acoustic model, the local 
peaks do not match with those of the compressible model. 

The sharp peak of the pressure amplitude near a frequency ratio of / - 0.0 is 
inherent in all the models (compressible, incompressible, and acoustic), shown in 
Figure 25. The phenomenon causing this resonant excitation is not yet understood and 
cannot be explained. However, the fact that this resonance occurs in all three models 
explains that neither the compressibility or the acoustics effects of the models affects 
this resonant frequency. Also, the fact that the frequency at which resonance occurs 
is near / = 0.0, indicates that this peak could be the response to a free vibration of 


P i 


=0 


(60) 
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the fluid inside the leakage path. 

Figure 26 shows the first-order pressure response for multiple cases of n. The 
frequency responses due to different n show no quantitative trend with respect to n. 
Unlike the compressible model, where the pressure response showed a definite, 
noticeable response to different n, the acoustic model shows little quantitative and 
qualitative response to different n value. 
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Figure 25 - Pressure response of incompressible, compressible, and acoustic model 
for impeller discharge pressure excitation for n = 1 for u m (0) = 0.7 
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CHAPTER VIII 

SUMMARY AND CONCLUSIONS 

A modified approach to Childs' previous work (1989,1992) on fluid-structure 
interaction forces in the leakage path between an impeller shroud and its housing is 
presented in this thesis. Three governing equations consisting of a continuity, path- 
momentum, and circumferential-momentum equations were developed to describe the 
leakage path inside a pump impeller. Radial displacement perturbations were used to 
solve for radial and circumferential force coefficients. In addition, impeller-discharge 
pressure disturbances were used to obtain pressure oscillation responses due to 
precessing impeller pressure wave pattern. Childs' model was modified from an 
incompressible model to a compressible barotropic-fluid model (the density of the 
working fluid is a function of the pressure and a constant temperature only). Results 
obtained from this model yielded interaction forces for radial and circumferential force 
coefficients. Radial and circumferential forces define reaction forces within the 
impeller leakage path. 

An acoustic model for the same leakage path was also developed. The 
convective, Coriolis, and centrifugal acceleration terms are removed from the 
compressible model to obtain the acoustics model. The compressible model is 
compared with the incompressible model and the acoustic model. A solution due to 
impeller discharge pressure disturbances model was also developed for the 
compressible and acoustic models. The results from these modifications are used to 
determine what effects additional perturbation terms in the compressible model have 
on the acoustic model. 

The results show that the additional fluid mechanics terms in the compressible 
model do cause resonances (peaks) in the force coefficient response curves. However, 
these peaks only occurred at high values of inlet circumferential velocity ratios. The 
peak pressure oscillation was shown to occur at the wearing ring seal. Introduction of 
impeller discharge disturbances with n diametral nodes showed that maximum peak 
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pressure oscillations occurred at nondimensional excitation frequencies of / = 6.4 and 
y _ 7 g f or this particular pump. Bolleter’s results suggest that for peak pressure 
oscillations to occur at the wearing ring seal, the nondimensional excitation frequency 
should be on the order of / = 2.182 for n = 11. The resonances found in this 
research do not match those predicted by Bolleter. At the predicted frequencies given 
by Bolleter, the compressible model shows an attenuation of the pressure oscillations 
at the seal exit. This does not provide a plausible explanation for the unusual behavior 
exhibited by the wearing ring seal. 

The compressibility of the fluid does not have a significant influence on the 
model at low values of nondimensional frequency. At high values of nondimensional 
frequency, the effects of compressibility become more significant. For the acoustic 
analysis, the convective, Coriolis, and centrifugal acceleration terms do affect the 
results to a limited extent for a precession excitation and a large extent for a pressure 
excitation when the fluid operates at relatively high mach numbers. 
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APPENDIX A 
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APPENDIX B 


Coefficient definitions for [B] of Eq. (50) 
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APPENDIX C 


Coefficient definitions for [C] of Eq. (57) 
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APPENDIX D 


Coefficient definitions for [D] of Eq. (60) 
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